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Abstract 

Branch  point  tolerant  phase  reconstructors  can  vastly  improve  adaptive  optic  system  performance  in 
extended  atmospheric  turbulence.  This  thesis  explores  the  performance  bounds  of  two  such  reconstructors 
Goldsteins  algorithm  and  hidden  phase.  A  least  squares  reconstructor  is  implemented  for  comparison. 
System  performance  is  presented  for  various  scenarios,  including  correction  time-delays,  wave-front 
sensor  noise,  and  extended  beacons.  These  scenarios  are  of  interest  for  laser  communication  and  directed 
energy  systems  such  as  Airborne  Laser.  Performance  bounds  are  obtained  through  wave-optics  simulation. 
The  extended  beacon  propagation  geometry  approximates  the  USAF  AFRL-DE  North  Oscura  Peak  range. 
Results  show  that  branch  point  tolerant  reconstructors  outperform  least  squares  for  equal  correction 
time-delays.  These  reconstructors  can  be  made  somewhat  tolerant  to  wave-front  sensor  error.  For  the  case 
of  an  incoherent  extended  beacon,  branch  point  information  is  lost  and  the  branch  point  algorithms 
perform  on  par  with  least  squares.  A  coherent  extended  beacon  preserves  branch  point  information,  but 
also  induces  branch  point  errors  due  to  coherent  speckle.  Still,  the  branch  point  reconstructors  tend  to 
maintain  a  1-2  order  of  magnitude  performance  advantage  over  least  squares  in  strong  turbulence.  While 
implementation  challenges  remain,  this  thesis  demonstrates  the  potential  of  branch  point  tolerant  phase 
reconstructors  on  laser  communication  and  weapons  systems. 
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Abstract 

Branch  point  tolerant  phase  reeonstructors  ean  vastly  improve  adaptive  optie  system 
performanee  in  extended  atmospherie  turbulenee.  This  thesis  explores  the  performanee 
bounds  of  two  such  reconstructors  -  Goldstein’s  algorithm  and  hidden  phase.  A  least 
squares  reeonstruetor  is  implemented  for  eomparison.  System  performanee  is  presented 
for  various  scenarios,  ineluding  correction  time-delays,  wave-front  sensor  noise,  and 
extended  beaeons.  These  seenarios  are  of  interest  for  laser  eommunieation  and  direeted 
energy  systems  sueh  as  Airborne  Easer.  Performanee  bounds  are  obtained  through  wave- 
opties  simulation.  The  extended  beaeon  propagation  geometry  approximates  the  EISAE 
AERE-DE  North  Oscura  Peak  range.  Results  show  that  braneh  point  tolerant 
reeonstruetors  outperform  least  squares  for  equal  eorreetion  time-delays.  These 
reeonstructors  can  be  made  somewhat  tolerant  to  wave-front  sensor  error.  Eor  the  ease  of 
an  ineoherent  extended  beaeon,  braneh  point  information  is  lost  and  the  braneh  point 
algorithms  perform  on  par  with  least  squares.  A  eoherent  extended  beaeon  preserves 
branch  point  information,  but  also  induces  braneh  point  errors  due  to  eoherent  speckle. 
Still,  the  braneh  point  reeonstruetors  tend  to  maintain  a  1-2  order  of  magnitude 
performanee  advantage  over  least  squares  in  strong  turbulenee.  While  implementation 
ehallenges  remain,  this  thesis  demonstrates  the  potential  of  braneh  point  tolerant  phase 
reeonstruetors  on  laser  eommunieation  and  weapons  systems. 
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DISTRIBUTED  BEACON  REQUIREMENTS  EOR  BRANCH  POINT  TOEERANT 
EASER  BEAM  COMPENSATION  IN  EXTENDED  ATMOSPHERIC  TURBUEENCE 

I.  Introduction 

Atmospheric  turbulence  is  the  primary  performance  limitation  for  earth  based 
optical  systems  that  operate  beyond  short  ranges.  Such  systems  include  astronomical 
telescopes,  optical  communication  systems,  and  directed  energy  weapons.  Through  the 
application  of  adaptive  optics  (AO),  the  effects  of  turbulence  can  be  reduced  and 
performance  enhanced.  The  degree  of  enhancement  depends  on  the  particular  system  and 
the  levels  of  turbulence  encountered. 

The  Airborne  Easer  (ABE)  represents  the  “state  of  the  art”  in  earth-based  directed 
energy  systems.  It  is  designed  to  engage  and  destroy  ballistic  missiles  at  ranges  of 
“hundreds  of  kilometers”  [1]  using  a  chemical  oxygen  iodine  laser  (COIL).  Adaptive 
optics  and  a  mature  COIL  design  are  the  foundations  for  ABL’s  shoot-to-kill  capability. 

Lour  lasers  and  a  deformable  mirror  comprise  the  major  optical  elements  of  the 
ABE.  Coarse  target  tracking  is  established  using  an  active  ranger.  The  target  is  then 
illuminated  with  the  track  illuminator  laser,  which  is  used  to  establish  the  aim  point  for 
the  COIL.  The  beacon  illuminator  laser  is  then  used  to  create  a  beacon  on  the  target  to 
measure  atmospheric  turbulence.  Ideally,  the  beacon  should  be  a  focused  spot  smaller 
than  the  diffraction  limit  of  the  ABE  system.  A  correction  is  calculated  and  applied  to 


the  deformable  mirror  from  which  the  COIL  kill  laser  is  propagated  to  the  target  and  held 
until  the  target  is  destroyed. 

The  adaptive  optics  portion  of  ABL  is  a  subset  of  the  major  optical  components, 
namely  the  beacon  illuminator  system  and  deformable  mirror.  This  system  measures  and 
compensates  for  atmospheric  turbulence.  The  process  of  measuring  turbulence  and 
creating  deformable  mirror  corrections  is  called  phase  reconstruction.  Successful  phase 
reconstruction  allows  the  COIL  kill  laser  to  maintain  focus  and  thus  energy  delivery  over 
long  distances. 

The  engagement  range  for  ABL  is  thus  directly  related  to  the  quality  of  phase 
reconstruction.  Atmospheric  turbulence  strength  increases  with  engagement  range. 

Above  some  turbulence  strength,  a  phenomena  known  as  branch  points  [2]  occurs  in  the 
beacon  field  which  drastically  limits  conventional  phase  reconstruction  techniques. 

Strong  turbulence  also  affects  the  beacon  -  enlarging  its  extent  beyond  the  diffraction 
limit  of  ABL.  The  resulting  “extended  beacon”  induces  measurement  errors  which 
increase  with  beacon  extent.  The  extended  beacon  may  also  become  scintillated  or 
broken  up  into  intensity  hot  spots  at  the  target. 

In  spite  of  the  branch  point  and  extended  beacon  effects,  ABL  is  projected  to  meet 
its  performance  goals  using  conventional  adaptive  optics  [1].  This  thesis  establishes 
performance  bounds  on  the  use  of  branch  point  tolerant  phase  reconstruction  in  the  ABL 
context  to  provide  guidance  on  future  ABL  improvements.  Since  extended  beacons  are 
expected  in  strong,  extended  turbulence,  this  thesis  also  considers  their  impact  on  branch 
point  reconstruction  performance.  Phase  reconstruction  is  primarily  a  software  function. 
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thus  these  improvements  eould  possibly  be  applied  without  ehanging  the  existing 
hardware. 

The  use  of  braneh  point  reeonstruetors  [3-5]  offers  two  eonsiderable  advantages 
over  the  current  conventional  reconstruction  [1],  First,  the  engagement  range  of  ABL  can 
be  extended  providing  a  larger  radius  of  theater  defense  and  allowing  more  flexibility  in 
protecting  the  ABL  platform  from  hostile  threats.  For  many  engagements  within  the 
current  envelope,  branch  point  reconstruction  will  place  more  energy  on  target,  reducing 
the  required  COIL  dwell  times  to  achieve  kill.  Since  COIL  requires  chemical  fuel,  a 
reduction  in  dwell  time  translates  into  more  shots  that  can  be  taken  before  the  ABL  must 
leave  station  for  replenishment.  For  follow-on  systems  to  the  ABL,  the  use  of  branch 
point  tolerant  reconstruction  may  allow  for  lower  power  requirements  in  the  kill  laser, 
probably  reducing  weight,  cost,  and  complexity. 

Atmospheric  turbulence  models  are  generally  characterized  as  either  “strong”  or 
“weak”  [6].  Weak  turbulence  conditions  allow  for  many  assumptions  that  simplify  the 
implementation  of  practical  AO  systems.  One  predominant  assumption  in  weak 
turbulence  is  the  absence  of  scintillation  or  amplitude  effects  on  the  optical  field.  This 
assumption  is  known  as  the  “near-field”  or  “phase-only”  turbulence  condition  since  the 
turbulence  can  be  modeled  as  a  thin  layer  of  phase  perturbation  close  to  the  receiving 
aperture.  Phase-only  conditions  do  not  affect  the  phase  continuity  of  incident  optical 
fields.  The  presence  of  continuous  phase  at  the  receiver  is  a  requirement  for  accurate 
least  squares  phase  reconstruction,  the  most  common  and  mature  optical  phase 
reconstruction  technique  [7].  Most  ground-based  astronomical  imaging  applications  can 
be  designed  assuming  phase-only  effects  [8]. 
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Terrestrial  eommunieation  and  direeted  energy  systems  often  must  operate  in  the 
presenee  of  phase  and  amplitude  effeets  [9].  As  first  discussed  by  Fried  [3],  amplitude 
effects  create  “branch  points”  or  discontinuities  in  the  optical  phase.  Least  squares 
techniques  average  out  these  discontinuities  as  noise  and  the  resulting  phase 
reconstruction  is  overly  smooth  [7].  The  performance  of  least  squares  AO  systems 
quickly  degrades  with  increasing  turbulence  strength  as  the  phase-only  assumptions  are 
violated  [10]. 

Various  techniques  have  been  proposed  to  develop  branch  point  tolerant  AO 
systems.  One  of  the  first  was  Fried’s  development  of  an  expression  for  the  “hidden 
phase”  or  phase  neglected  by  least  squares  algorithms,  which  could  then  be  summed  with 
a  least  squares  output  and  used  for  correction  [3].  Parallel  to  this  approach  is  the  method 
of  slope  discontinuity  -  wherein  the  phase  is  expressed  as  the  sum  of  a  vector  curl  and 
scalar  divergence  representing  the  least  squares  and  hidden  phase  components  [5].  The 
slope  discontinuity,  though  analogous  to  hidden  phase,  does  not  require  that  actual  branch 
point  locations  be  determined.  It  is  generally  believed  that  determining  branch  point 
locations  using  practical  wave  front  sensors  is  not  realistic  at  present  [11]. 

While  work  is  under  way  developing  realizable  branch  point  tolerant 
reconstructors,  it  is  still  important  to  quantify  the  upper  bounds  of  achievable 
performance.  The  use  of  Goldstein’s  algorithm  [7]  by  Roggemann  was  one  of  the  first 
such  studies  [4].  Roggemann  found  that  Goldstein’s  algorithm  provided  good  upper 
bounds  for  noise-free  measurements  and  zero  time-delay  corrections  using  a  finite 
resolution  deformable  mirror.  Goldstein’s  algorithm  is  a  branch  point  tolerant  phase 
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reconstructor  developed  for  use  in  synthetic  aperture  radar,  but  is  also  applicable  to 
optical  phase  reconstruction. 

Branch  point  reconstructor  studies  to  date  [3-5]  have  assumed  point  source 
beacons.  This  assumption  is  valid  whenever  the  finite  beacon  extent  is  less  than  the 
diffraction-limited  resolution  of  the  optical  system.  Unfortunately,  in  the  Airborne  Laser 
and  other  cases,  the  beacon  extent  often  exceeds  the  diffraction-limited  resolution  and 
induces  errors  in  phase  measurement.  Such  a  beacon  is  defined  as  an  “on-axis  circular 
extended  beacon  [12].”  This  beacon  is  a  specific  case  of  a  distributed  beacon.  In 
general,  distributed  beacons  violate  ideal  beacon  assumptions  by  some  combination  of 
off-axis  location  and  finite  extent. 

This  thesis  establishes  the  constraints  on  extended  beacon  radius  in  order  to 
implement  branch  point  tolerant  reconstructors  successfully.  This  is  done  through  wave- 
optics  simulation.  Simulations  are  conducted  using  the  Air  Force  Research  Labs  North 
Oscura  Peak  (NOP)  geometry.  NOP  is  a  specialized  adaptive  optics  range  in  the  White 
Sands  Missile  range,  which  has  been  used  extensively  for  Airborne  Laser  development. 
While  the  NOP  geometry  is  not  an  exact  match  to  ABL,  results  from  NOP  are  directly 
scalable  and  thus  applicable  to  ABL  scenarios. 

Two  forms  of  extended  beacons  are  considered  in  simulation  space-  a  coherent 
extended  beacon  and  an  incoherent  extended  beacon.  These  represent  the  limiting 
extremes  of  beacon  coherence  in  an  adaptive  optics  system.  Performance  bounds  for 
various  reconstructor  algorithms  are  established  for  various  beacon  sizes  and  turbulence 
strengths.  Point  source  performance  is  also  presented  to  establish  a  basis  for  comparison. 
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This  thesis  is  organized  as  follows.  Chapter  2  provides  a  review  of  published 
literature  on  phase  reeonstruetion  in  general  and  braneh  points  specifieally.  A  discussion 
of  atmospheric  turbulence  modeling  follows  in  Chapter  3.  This  chapter  also  establishes 
the  theoretical  concepts  necessary  for  implementation  of  reconstruction  techniques  in 
practice  or  simulation.  Chapter  4  provides  specific  details  on  the  simulations  conducted 
for  this  thesis.  These  details  include  considerations  for  choosing  modeling  parameters 
and  the  construction  of  extended  beacon  models.  Simulation  results  are  provided  in 
Chapter  5,  along  with  discussion  and  comments.  Conclusions  are  presented  in  Chapter  6 
along  with  recommendations  for  future  research. 
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II.  Literature  Review 


This  chapter  presents  a  literature  review  on  phase  reeonstruction  and  distributed 
beacon  research  relevant  to  this  thesis.  Three  phase  reeonstruetors  are  discussed  -  least 
squares,  hidden  phase,  and  Goldstein’s  algorithm.  References  are  also  provided  on  past 
work  with  distributed  beaeons. 

2.1  Least  Squares  Reconstruction 

Adaptive  optics  techniques  to  correct  for  near-field  or  phase-only  turbulenee  are 
well  developed  [8].  They  have  seen  wide  applieation  in  astronomieal  imaging. 

Generally,  the  least  squares  reconstruetor  is  used  for  phase  reeonstruction  and  is 
developed  from  the  following  equation  [4]; 

A\j7  =  P\j7  (1) 

where  A\j7  represents  a  vector  of  wrapped  x  and  y  phase  differences  taken  in  a  two- 
dimensional  plane  like  those  from  a  wave-front  sensor,  tj7  is  a  matrix  of  input  phase 
values,  and  P  is  a  mapping  from  \j7  to  A\j7  .  For  a  given  input  phase  of  dimension  MxN, 
we  will  have  Mx(N-l)  phase  differenees  in  x  and  (M-l)xN  phase  differenees  in  y.  If  we 
had  only  MxN  phase  differences,  the  straight  forward  solution  to  Equation  (1)  would  be 
tj;  =  P'^Atj/ ,  but  A\j/  is  comprised  of  nearly  (MN)^  differences  so  the  order  of  P'^Atj/ 
exceeds  that  of  tj/  and  a  singular  solution  is  not  possible.  As  is  often  the  case,  we  opt  for 
an  L  minimized  solution  for  sueh  an  over-eonstrained  problem.  This  can  be  given  as 

vi/„  =  (P"P)-^P"Atj/  (2) 
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where  represents  the  least  squares  reeonstruetion  of  the  wrapped  phase  differenees 
Atj/  .  The  form  of  P  is  important  in  determining  how  best  to  solve  Equation  (2). 
Consider  the  simple  example  of  a  3x3  grid  of  input  phase  measurements 


¥n  ¥21  ¥31 

¥12  ¥22  ¥32 

¥13  ¥23  ¥33 


(3) 


There  are  6  possible  Ax  phase  measurements  and  6  possible  Ay  phase  measurements 
between  adjaeent  inputs.  Let  A¥"  be  a  12-element  vector  containing  the  6  Ax  slopes 
followed  by  the  6  Ay  slopes.  We  wish  to  map  the  3x3  grid  above  to  A¥^  via  a  matrix  P. 
We  cannot  go  directly  from  a  3x3  matrix  to  a  1x12  vector,  so  first,  we  reshape  the  input 

matrix  into  a  9x1  vector  ¥^  =  [ ¥11  ¥21  ¥31  ¥12  ¥22  ¥32  ¥13  ¥23  ¥33]^-  Now, 
we  can  define  a  matrix  P  that  is  12x9  to  accomplish  our  mapping.  The  resulting  matrix  is 


-1  1  0  0  0  0  0  0 

0  -1  1  0  0  0  0  0 

0  0  0  -1  1  0  0  0 

0  0  0  0  -1  1  0  0 

0  0  0  0  0  0  -1  1 

0000000  -1 
-1  0  0  1  0  0  0  0 

0  -1  0  0  1  0  0  0 

0  0  -1  0  0  1  0  0 

0  0  0  -1  0  0  1  0 

0  0  0  0  -1  0  0  1 

00000  -1  00 


0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

1 


(4) 


which  is  composed  entirely  of  I’s,  -I’s,  and  O’s  and  for  large  input,  grids  will  be  sparse. 
Efficient  solutions  of  Equation  (2)  have  been  developed  using  iterative,  sparse  matrix, 
and  transform  techniques  [7].  The  performance  of  least  squares  systems  in  weak 
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turbulence  can  be  very  good,  but  the  presence  of  branch  points  drastically  reduces  their 
effectiveness,  as  shown  by  Primmermann  in  1995  [10]. 

2,2  Branch  Points 

Atmospheric  turbulence  gives  rise  to  random  perturbations  in  the  amplitude  and 
phase  of  an  optical  field.  In  strong  turbulence,  amplitude  effects  can  cause  the  optical 
field  to  have  localized  points  of  zero  amplitude.  At  such  points,  the  phase  is  undefined 
and  this  gives  rise  to  a  discontinuous  phase  function.  These  points  are  called  branch 
points  [2]  and  they  cause  phase  reconstruction  errors  when  using  least  squares 
reconstruction.  The  idea  of  branch  points  and  discontinuous  phase  was  first  discussed  in 
an  optics  setting  by  Fried  in  1992  [2]. 

A  1998  paper  by  V.  Aksenov  et.  al.  provided  one  of  the  first  mathematical  models 
for  branch  points  in  optical  fields  [13].  This  paper  noted  that  a  discontinuous  phase 
exhibited  both  potential  and  vortex  features.  Using  Akesenov’s  paper  as  a  starting  point. 
Fried  developed  a  two-component  phase  model.  This  model  expressed  the  phase  as  the 
sum  of  a  scalar  gradient  and  vector  curl.  We  will  recreate  Fried’s  development  in  [3] 
beginning  with  his  Equation  (7), 

±271  if  branch  point  enclosed 

(5) 

0  if  no  branch  point  enclosed 
where  g(r(^))  is  the  wrapped  phase  gradient  of  the  input  phase  at  a  point  r(^)  defined  on 
the  contour  C  and  t(^)  defines  a  unit  vector  tangent  to  the  contour  at  ^  times  some 
length  coefficient.  Fried  then  invokes  Stake’s  theorem  to  express  Equation  (5)  as 
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■  Vxg(r) 


±271  if  branch  point  enclosed 
0  if  no  braneh  point  enelosed 


(6) 


where  \  is  a  unit  vector  in  the  z  direetion.  He  then  eonsiders  the  limiting  ease  of  a  small 
eontour  enelosing  a  single  braneh  point.  The  integrand  then  reduces  to  a  Dirae  delta 
function  centered  at  the  branch  point  location  r^p  .  The  result  from  Equation  (6)  is 

l,-Vxg(r)  =  ±27t5(r-r,p).  (7) 

We  know  that  the  eurl  of  a  gradient  is  identieally  zero,  so  g(r)  must  represent  more  than 
just  the  gradient  of  a  scalar  phase  function.  Since  g(r)  is  a  vector  function,  it  can  be 
expressed  as  a  sum  of  a  sealar  gradient  and  the  eurl  of  a  veetor  funetion  as  follows: 

g(r)  =  V5(r)  +  VxH(r).  (8) 

We  reeall  that  the  least  squares  solution  provided  a  path-independent  mapping 
from  sealar  phase  differenees  (analogous  to  gradients)  to  unwrapped  input  phase.  These 
sealar  phase  differenees  are  represented  by  V5(r)  in  Equation  (8).  When  braneh  points 
are  present,  H(r)  is  non-zero  and  the  least  squares  reeonstruction  is  missing  information. 
Intuitively,  the  braneh  points  eause  the  reeonstruetion  to  have  some  path  dependenee  -  a 
path  which  crosses  a  line  between  branch  points  gets  some  ±27i  jump.  These  lines  are 
ealled  braneh  cuts.  Unlike  typieal  measurement  noise,  the  phase  jumps  are  not  zero- 
mean  and,  when  averaged,  will  cause  reeonstruetion  errors.  Eurther,  the  phase  error  is 
not  restrieted  to  the  region  near  the  braneh  cut;  a  phase  jump  affeets  every  point  in  the 
grid  sinee  least  squares  eonsiders  all  possible  paths.  Two  solutions  logieally  present 
themselves  -  develop  some  solution  for  VxH(r)  whieh  ean  eorreet  the  least  squares 
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reconstruction  or  locate  and  avoid  branch  cuts  altogether  and  perform  a  path  dependent 
reconstruction. 

2,3  Hidden  Phase  Reconstructor 

Fried  chose  the  former  solution  and  formulated  an  expression  for  “hidden  phase.” 
Hidden  phase  is  defined  as  the  phase  error  represented  in  neglecting  VxH(r)  during 
least  squares  reconstruction.  Returning  to  Equation  (7)  and  (8),  we  remember  that  the 
curl  of  the  scalar  divergence  is  identically  zero.  Substituting  Equation  (8)  into 
Equation  (7), 

l^-VxVxH(r)  =  -V"/2(r)  (9) 

where  /i(r)  is  the  z  component  of  H(r)  [3].  Relating  Equation  (9)  to  Equation  (7), 

V^/i(r)  =  +27i5(r-r,p)  (10) 

which  has  a  solution  given  by 

/i(r)  =  +log(|r-r^p|).  (11) 

Having  already  related  the  least  squares  solution  to  the  scalar  divergence  operator  in 
Equation  (8),  Eried  equates  the  divergence  of  the  hidden  phase  to  VxH(r)  and  after 
some  manipulation,  is  able  to  show  that 

^ud (r)  =  Im{±log [(x - +  / [y  - )]}  .  (12) 

Hidden  phase  is  then  calculated  by  evaluating  \|/^.^  (r)  for  all  r  for  each  branch  point 
location  {xbp,  ytp)  and  summing  the  result.  The  total  phase  is  then  expressed  as 

^=^LS+^hid  (13) 
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Since  the  hidden  phase  operates  in  parallel  with  a  least  squares  reconstruction  and 
requires  the  same  inputs,  implementation  is  straightforward  [3]. 

Fried’s  paper  was  soon  followed  by  Le  Bigot  and  Kibblewhite  who  proposed  a 
method  for  braneh  point  deteetion  based  on  an  iterative,  intensity-weighted  reeonstruetor 
[14],  Voitsekhovich  et.  al.  published  a  lengthy  paper  on  braneh  point  densities  [15].  This 
paper  presented  both  theoretical  and  experimental  branch  point  densities  based  on 
turbulenee  strength.  As  expeeted,  it  showed  that  braneh  points  inereased  as  a  funetion  of 
propagation  distance  and  scintillation  index. 

Le  Bigot  and  Wild  collaborated  on  a  1999  study  which  proposed  yet  another 
method  for  branch  point  detection  [11].  This  method  included  the  modeling  of 
measurement  noise.  The  same  year,  Arrasmith  published  the  first  results  from  a  hidden 
phase  reeonstruetor  implementation  [16].  He  obtained  promising  results  using  the  hidden 
phase  reeonstruetor  in  eoherent  image  restoration  and  did  eonsider  measurement  noise. 

2.4  Goldstein’s  Algorithm  Reeonstruetor 

The  first  use  of  Goldstein’s  reeonstruetion  algorithm  [7]  in  an  adaptive  optics 
settings  was  documented  by  Roggemann  and  Koivunen  in  2000  [4].  While  Goldstein’s 
algorithm  is  difficult  to  implement  in  praetice,  it  provides  exeellent  upper  bounds  on 
branch  point  tolerant  reeonstruetor  performance.  Goldstein’s  algorithm  uses  branch  point 
locations  to  place  braneh  euts  and  effeet  a  path-dependent  reeonstruetion  around  the  cuts. 
Roggemann’ s  study  was  one  of  the  first  to  examine  branch  point  reeonstruetor 
performance  extensively  in  an  extended  turbulenee  setting. 
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Goldstein’s  algorithm  is  explained  in  detail  by  [7],  but  it  is  important  to 
understand  the  mechanism  behind  branch  cut  creation.  A  key  element  in  placing  branch 
cuts  is  the  minimization  of  cut  length.  Cut  placement  begins  with  a  search  of  the 
wrapped  phase  for  branch  points.  When  a  branch  point  is  found,  the  algorithm  searches 
the  3x3  neighborhood  of  the  branch  point  for  another  branch  point.  If  another  branch 
point  is  found,  a  cut  is  placed  between  the  two.  Should  the  branch  points  be  a  positive¬ 
negative  pair,  they  are  designated  “balanced”.  If  the  branch  points  are  of  the  same  type, 
the  neighborhood  of  the  branch  point  is  then  searched  for  other  branch  points.  A 
branch  cut  is  considered  balanced  then  the  sum  of  positive  and  negative  branch  points  is 
equal.  While  a  cut  is  unbalanced,  the  algorithm  will  continue  to  shift  and  grow  the  search 
neighborhood  until  balance  is  achieved  or  the  image  edge  is  reached.  If  the  edge  is 
reached  first,  the  cut  will  terminate  at  the  border  and  the  cut  will  be  considered  balanced. 
Once  all  the  cuts  are  in  place,  the  phase  is  unwrapped  around  the  cuts  using  the  wrapped 
phase  inputs. 

Recognizing  the  limitations  in  implementing  Goldstein’s  algorithm  on  real 
systems,  Roggemann  and  Koivunen  followed  up  with  another  publication  which  outlined 
the  design  of  an  iterative  intensity-weighted  least  squares  reconstructor  which  approached 
Goldstein’s  performance  without  the  need  to  find  branch  point  locations  [17].  Another 
branch  point  tolerant  reconstructor  was  proposed  by  Tyler  [5].  He  formulated  the  phase 
as  a  summation  of  the  least  squares  solution  with  slope  discrepancy.  The  slope 
discrepancy  was  a  reformulation  of  Fried’s  hidden  phase  without  the  need  for  branch 
point  locations  and  is  somewhat  similar  to  Le  Bigot  and  Wild’s  ideas  in  1999.  These 
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reconstructors  are  mentioned  for  completeness  and  were  not  considered  as  part  of  this 
thesis. 

2,5  Extended  Beacons 

While  branch  point  reconstructors  have  been  shown  to  work  well  for  point  source 
beacons,  little  has  been  published  regarding  their  application  to  extended  beacons.  If 
branch  point  reconstructors  are  to  be  utilized  on  ABL  or  other  directed  energy  systems 
which  operate  in  strong  turbulence,  they  must  be  able  to  function  with  extended  beacons. 
The  theory  of  extended  beacons  has  its  roots  in  studies  conducted  for  artificial  guide  stars 
and  active  illumination  [18-20].  These  studies  examined  the  effects  of  various  forms  of 
anisoplanatism  (off-axis  beacons)  for  different  adaptive  optics  applications.  Stroud 
expanded  on  these  studies  to  include  an  analysis  of  the  effect  of  finite  extent  [12]  for  on 
on-axis  uniform  intensity  circular  beacon.  He  found  that  the  on-axis  extended  beacon 
had  the  least  performance  degradation  of  the  types  of  distributed  beacons  he  studied.  Rao 
also  found  that  the  on-axis  extended  beacon  was  the  most  favorable  distributed  beacon 
[21]. 


We  have  now  established  the  necessary  background  for  extended  beacons  and 
phase  reconstruction.  With  this  background,  we  are  ready  to  formulate  simulations  to 
measure  the  performance  of  the  phase  reconstructors  in  directed  energy  applications. 
Chapter  3  will  present  the  necessary  theory  to  conduct  simulations  with  atmospheric 
turbulence  modeling  using  the  least  squares,  hidden  phase,  and  Goldstein  reconstructors. 
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III.  Theory 


Theory  for  adaptive  optics  can  be  divided  into  two  broad  categories  -turbulence 
modeling,  and  AO  system  modeling.  In  the  context  of  this  research,  turbulence  modeling 
describes  the  effects  of  the  earth’s  atmosphere  on  laser  propagation  while  AO  system 
modeling  is  focused  on  phase  reconstruction  from  beacon  imagery.  Since  we  are  dealing 
with  simulations  and  not  real  hardware,  most  other  issues  in  the  AO  system  modeling 
(mirror  alignments,  platform  jitter,  heat  dissipation,  etc. . .)  can  be  neglected  as  outside  the 
scope  of  investigation. 

3.1  Atmospheric  Turbulence  Modeling 

Atmospheric  turbulence  -  both  the  weather  and  optical  varieties  -  is  caused  by 
uneven  solar  heating  of  the  Earth’s  surface.  The  temperature  dependence  of  the  index  of 
refraction  is  given  by 

n  =  1  +  77.6(1  +  7.52x10"' /l"")^xlO"®  (14) 

where  P  is  the  atmospheric  pressure  in  millibars,  T  the  temperature  in  degrees  Kelvin, 
and  X  the  wavelength  of  light  in  micrometers  [22].  Changes  in  n  are  dominated  by 
temperature  changes,  and  X  and  P  are  normally  treated  as  constants. 

Given  the  relationship  between  temperature  and  index  of  refraction,  the  next  step 
towards  an  atmospheric  model  is  to  determine  the  distribution  or  power  spectral  density 
of  these  index  variations  in  the  atmosphere.  Uneven  solar  heating  creates  large  zones  of 
similar  temperature.  These  zones  or  “turbulent  eddies”  [9]  are  broken  up  by  wind  and 
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convection  into  ever  smaller  regions  of  homogeneous  index  of  refraction.  The  power 
spectral  density  represents  the  relative  quantity  of  eddies  with  a  given  dimension 

=  2n ! ,  Ly  =  27[ I Ky,  and  1^=271 ! In  optics,  we  generally  assume  isotropic 
turbulence  which  reduces  the  dependence  of  to  L  =  2n !  k  where  /r  =  |^| . 

Based  on  Kolmogorov  theory  [23],  is  generally  described  by  three 

regions.  For  /rless  than  some  value  (large  scale  sizes),  the  power  spectral  density  is 
dominated  by  large-scale  geographic  and  meteorological  forces,  which  cannot  be 
predicted  by  turbulent  flow  theory.  The  dimension  L^=27l !  is  called  the  “outer 
scale”  and  marks  the  boundary  dimension  below  which  the  form  of  behaves 

according  to  turbulence  theory.  Lo  is  typically  on  the  order  of  1  to  100m,  depending  on 
geometry  and  atmospheric  conditions.  This  region  of  /r  >  is  known  as  the  inertial 

subrange  and  is  based  on  Kolmogorov’s  theory  of  turbulent  flows.  The  PSD  can  then  be 
expressed  as 

4)„(x-)  =  0.033C„V-"''  (15) 

where  C„  is  the  atmospheric  structure  constant  of  the  index  variations  and  is  a  general 
indication  of  turbulence  strength. 

Predictably,  there  is  also  an  upper  bound  on  K,  Km,  above  which  <I>„(/f)  drops 
rapidly  as  turbulent  eddies  dissipate  through  viscous  forces.  The  dimension  1^=271 ! 
or  inner  scale  describes  this  limit.  The  form  of  <I>„(/f')  for  the  Kolmogorov  spectrum 
given  in  Equation  (15)  is  thus  valid  over  k^<k  <K^,ox  as  more  commonly  written, 

I/Iq  «  1/4 . 
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Various  extensions  to  the  Kolmogorov  speetrum  have  been  made  to  model 
behavior  outside  the  inertial  subrange.  The  Tatarskii  speetrum  [24]  extends  into 


the  dissipation  range  by  adding  a  truneating  funetion  to  Equation  (15), 


<D„(x-)  =  0.033C„V-"''exp - ^  ,  x-»l/4. 


The  Tatarskii  speetrum  ean  be  further  enhaneed  by  extension  into  the  region  K<  ULo. 

The  most  eommon  model  of  this  type  being  the  von  Karman  speetrum.  While  the  von 
Karman  and  other  turbulenee  models  [25-27]  are  more  aeeurate  over  large  wave 
numbers,  the  Kolmogorov  and  Tatarskii  speetrums  are  eommonly  used  whenever 
possible  due  to  their  simplieity. 

In  praetieal  situations,  C„  often  varies  with  propagation  distanee.  Most 
fluetuations  are  due  to  slant  or  vertieal  paths  through  the  atmosphere  as  eneountered  in 
astronomy,  eommunieations,  and  direeted  energy  applieations.  For  these  eases,  the 
power  speetral  density  ean  be  rewritten  as 

<I>„(x-,z)  =  0.033C„"(z)x--"'^  (17) 

to  eapture  the  dependenee  of  C„  on  propagation  distanee.  Various  models  have  been 

2 

developed  based  on  atmospherie  measurements  for  C„  vs.  altitude  [28]. 


Atmospheric  turbulence  strength  is  generally  characterized  by  the  Rytov  variance 
or  log-amplitude  variance  of  the  incident  field.  Theoretically,  Rytov  variance  is 
calculated  using  the  following  equation  for  a  spherical  wave  propagation  through 
turbulence  [28], 

L  /  \5/6 

\dzCl{z)  -  {L-zf"  (18) 

0  v^7 
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where  A:  =  2;r/  A ,  X  is  the  wavelength,  L  is  the  propagation  distance,  andC„^  (z)  is  the 
atmospheric  structure  constant.  The  Rytov  variance  is  also  commonly  referred  to  as  the 
Rytov  number  or  just  Rytov.  For  the  case  where  Cl  (z)  is  constant  over  the  propagation 
path,  the  integral  simplifies  to 

(19) 


Since  simulations  are  discrete,  we  must  break  the  turbulence  path  into  segments 
or  “slabs”  of  turbulence.  We  generally  calculate  the  total  phase  perturbation  in  each  slab 
and  apply  it  to  the  propagating  field  via  thin  phase-screens  [8].  For  the  complex 
representation  of  U(x),  we  express  the  incident  phase  as 


=  tan 


'lm[17(x)]~ 

Re[f/(x)] 


(20) 


Phase-only  turbulence  is  then  applied  via 

(x)  =  (x)  exp  [jy/^  (x)] , 


(21) 


where  U ^  (x)  is  the  field  incident  on  a  phase  screen  and  xj/n  (x)  is  the  phase  perturbation 
applied  by  the  screen.  Screens  are  placed  a  sufficient  distance  apart  such  that  they  are 
statistically  independent  of  each  other.  This  simplifies  implementation  by  removing  the 
need  to  account  for  correlations  between  phase  screens.  The  Rytov  variance  of  N  screens 
is  calculated  by  discretizing  Equation  (18), 

(22) 

«=i  V  y 

where  is  the  distance  between  each  screen.  We  can  also  think  of  as  the  slab 
thickness  for  each  of  our  discretized  turbulence  segments. 
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When  the  Rytov  variance  exceeds  about  0.05,  both  amplitude  (scintillation)  and 
phase  effects  are  generally  present  in  the  propagation  path  [4].  Scintillation  indicates 
branch  points  in  the  phase.  It  is  appropriate  to  examine  some  samples  of  point  source 
propagations  taken  over  various  turbulence  strengths.  Figure  1  shows  the  intensity  and 
phase  for  a  50  km  point  source  propagation  in  a  vacuum.  We  see  that  the  intensity  yields 
a  diffraction  limited  Gaussian  spot  and  that  the  phase  changes  slowly  over  the  aperture 
with  a  well-defined  structure. 


Figure  1:  Intensity  and  Phase  for  50  km  Point  Source  Propagation  in  Vacuum 
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In  Figure  2,  we  begin  to  see  the  onset  of  scintillation  in  moderate  turbulence.  The 
intensity  has  been  broken  up  into  3  definite  hot  spots  and  some  energy  has  been  spread 
across  the  aperture.  These  hot  spots  are  scintillation.  The  phase  map  also  reveals  that  the 
phase  is  changing  more  rapidly  across  the  aperture  and  the  structure  is  becoming  more 
random.  Finally,  in  Figure  3,  we  see  the  effect  of  strong  turbulence  on  the  intensity  and 
phase.  The  intensity  has  been  further  spread  throughout  the  aperture,  but  the  change  in 
phase  is  even  more  dramatic.  In  this  case,  the  phase  structure  is  lost  and  there  are  rapid 
fluctuations  in  value  throughout. 


Figure  3;  Intensity  and  Phase  for  Point  Source  Propagation  in  Strong  Turbulence 


The  scintillated  fields  in  Figure  2  and  Figure  3  contain  branch  points.  These 
discontinuities  cause  errors  in  least  squares  reconstructions.  In  a  sampled  phase  field, 
branch  points  are  located  by  calculating  the  residuals  over  each  2x2  pixel  block  via  [3]: 

X ,,,  (h  j)  =  -A^  (/,  j)  -  A ”  (/  + 1,  j)  +  (i,  7  + 1)  +  A ’■  (/,  j)  (23) 

where  and  represent  the  phase  slopes  calculated  using  the  following 

equations  when  the  incident  phase,  y/j  (/,  j) ,  is  known; 
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A"  (/,  j)  =  arg(exp  (/  + 1,  j)  -  y/,  {i,  7)] } ) 


(24) 


(/,  j)  =  arg(exp  {i,  7  + 1)  -  (/,  7)] } )  (25) 

A  summation  in  Equation  (23)  equal  to  2%  or  -2%  indicates  the  presence  of  a 
positive  or  negative  branch  point.  In  the  sampled  case,  the  branch  point  location  can  only 
be  localized  to  within  a  2x2  grid  and  the  normal  practice  is  to  assign  its  location  to  the 
upper  left  pixel.  Clearly,  finite  sampling  can  have  a  large  impact  on  localizing  branch 
points  and  is  a  major  concern  for  the  practical  application  of  the  Goldstein  and  hidden 
phase  algorithms  [11]. 

3.2  Parameterization  of  Turbulence 

While  Rytov  variance  ( )  is  the  most  common  variable  for  atmospheric 

parameter  studies,  it  is  certainly  not  the  only  parameter  of  interest.  Two  other  common 
parameters  are  and  60.  The  parameter  ro  is  known  as  the  atmospheric  coherence 
diameter  or  Fried  Parameter  [9].  It  is  often  used  to  estimate  the  resolution  of  optical 
systems  in  turbulence.  6^  is  the  isoplanatic  angle  which  defines  an  angular  radius  from  a 
beacon  within  which  the  turbulence  effects  are  essentially  constant  [6].  The  general 
expression  for  G^  was  given  previously  in  Equation  (18).  The  expressions  for  Tq  and  9o 
are  as  follows: 
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0 

and  assume  spherical  wave  propagation  with  Kolmogorov  turbulence.  Discretized  for 
simulation  via  phase  screens,  these  equations  become 
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(28) 


<j|=0.563t”‘r'‘C„£A,Jf., 


n=\ 


f  2.9144  D 


(29) 


and 


(ej  =2.9144  eC,Y/3J„. 


(30) 


In  each  equation,  Co  is  a  nominal  turbulence  weight,  N  is  the  total  number  of 
phase  screens,  is  the  specific  turbulence  weight  for  the  n*  screen  and  T^,R„,  and 

represent  the  integrals  over  each  distance  interval  between  screen.  For  example,  in 
Equation  (28)  we  let 
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(31) 


where  r„  is  the  position  of  the  n  phase  screen.  Normally,  is  varied  by  changing  Co, 
which  also  causes  Tq  and  Og  to  vary.  This  is  the  important  result  to  consider; 
individually,  (7^ ,  ,  and  Og  cannot  uniquely  define  the  propagation  path  if  there  are  more 

than  3  phase  screens. 

If  we  desire  to  test  dependence  on  a  subset  of  <7^ ,  ,  or  6g,  we  must  do  a  bit  more 

work  in  determining  the  phase  screen  weights.  Since  Equations  (28)-(30)  are  linear  in 
,  we  can  represent  them  in  matrix  fashion: 

\<jl\=CgAP  with  A  \xN  (32) 


X 

, -5/3 


with  B  IxN 


(33) 
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=  c^cp  with  C  3xN  (34) 

and  let =  [y^j  ...  Pn^  ■  We  note  that  A,  B,  and  C  are  matriees  of  rank  3  and  null 

spaee  of  dimension  (N-3).  If  we  wish  to  eonduet  a  parameter  study  on  6^  while 
eonstraining  and  r^,  this  is  readily  aeeomplished  by  finding  null-spaee  veetors  in 
Equation  (33)  and  adding  them  to  P  in  Equation  (34).  Similarly,  with  slight 
modifieations,  the  same  teehnique  ean  be  used  to  eonstrain  any  desired  eombination  of 
the  three  variables  for  parametrie  studies. 

Parameterization  studies  eondueted  as  shown  above  are  best  implemented  when 
the  range  of  parameter  values  is  limited.  Physieal  turbulenee  sereens  eannot  have 
negative  strength,  whieh  limits  the  ehoiee  of  eoeffieients  for  the  null  spaee  veetors.  The 
work  in  this  thesis  is  eondueted  over  a  large  range  of  parameters.  To  avoid  long  searehes 
for  appropriate  null  spaee  veetors,  turbulenee  strengths  are  set  by  varying  Q  whieh 

results  in  the  simultaneous  variation  of  ,  and  O^. 

3.3  Reeonstruetor  Algorithms 

To  implement  adaptive  opties,  we  must  measure  the  effeets  of  the  atmosphere  and 
determine  an  appropriate  eorreetion.  Generally,  sueh  eorreetions  are  applied  via  a 
deformable  mirror.  Deforming  portions  of  the  mirror  will  induee  optieal  path  delays  on 
ineident  light  and  thus  ehange  the  phase  of  the  ineident  field.  As  just  seen  in  Eigure  1,  a 
diffraetion-limited  spot  produees  a  uniform  phase.  Eor  an  ideal  eorreetion  in  turbulenee, 
we  wish  to  apply  a  phase  eorreetion  sueh  that  the  final  phase  is  also  uniform.  Under  the 
assumption  that  the  beaeon  phase  was  uniform,  the  measured  phase  at  the  reeeiver 
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represents  the  atmospherie  distortion,  whieh  we  ean  eonjugate  and  apply  as  our 
eorreetion. 

To  do  so,  we  must  first  measure  the  ineident  wrapped  phase  and  apply  some  type 
of  phase  retrieval  algorithm.  Phase  eannot  be  measured  direetly  in  physieal  systems; 
rather,  eommon  wave-front  sensors  measure  the  phase  slope  or  tilt  over  a  finite  number 
of  subapertures  within  the  ineident  field.  Praetieal  optieal  reeonstruetor  algorithms  must 
be  formulated  with  the  use  of  phase  slopes  as  the  primary  input.  The  phase  slopes  are 
also  ealled  phase  differenees. 

The  least  squares  reeonstruetor  is  the  most  eommon  and  mature  in  adaptive 
opties  [8].  ft  is  readily  formulated  to  aeeept  phase  slopes  as  an  input  and  output  an 
estimate  of  the  unwrapped  phase.  The  least  squares  teehnique  unwraps  phase  aeross 
every  possible  path  whieh  tends  to  average  out  zero-mean  measurement  noises  [7]. 
Unfortunately,  it  also  averages  out  braneh  outs  sinoe  it  assumes  oontinuous  phase. 

Fried  developed  an  expression  for  the  “hidden  phase,”  Equation  (12),  or  phase 
error  oaused  by  applying  least  squares  reoonstruotion  when  braneh  points  are  present  [3]. 
The  hidden  phase  algorithm  oaloulates  the  least  squares  error  from  the  braneh  point 
looations  disoovered  through  Equation  (23).  Generalizing  Equation  (12)  for  K  braneh 
points,  hidden  phase  is  given  by 


log 


W{x-x^)  +  j{y-y,) 


k=\ 


^{x-x[)  +  j{y-y[) 

k=\ 


(35) 


where  Kp  is  the  number  of  positive  braneh  points,  Kn  is  the  number  of  negative  braneh 
points,  X  andy  are  the  pixel  looations,  ix^,y^)  is  the  looation  of  the  positive  braneh 
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point,  and  {x[,y[)  the  location  of  the  negative  branch  point.  In  practice,  Equation  (35) 
can  cause  a  numerical  overflow  when  operating  on  large  grids  in  strong  turbulence.  To 
avoid  this  problem,  an  alternate  expression 


¥/„y(r)  =  Im 


log 


n 

j^(Kp-Kn)  k=\ 


{x-x,)  +  j{y-y^) 
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{x-x[)  +  j{y-y[) 


It 

k=\ 


N 


(36) 


can  be  used,  where  N  is  the  grid  width  in  pixels  (assuming  a  square  grid).  Equation  (36) 
eliminates  overflow  for  the  range  of  turbulence  strengths  and  grid  sizes  simulated  in  this 
thesis.  The  final  phase  used  for  reconstruction  then  becomes 

¥(r)  =  ¥/„.(r)  +  ¥is(r)  (37) 

where  rir^jlr)  is  the  reconstructed  phase  from  a  least  squares  algorithm. 

Another  branch  point  tolerant  phase  reconstructor  is  Goldstein’s  algorithm.  This 
algorithm  also  calculates  the  location  of  branch  points  through  Equation  (23).  The 
branch  points  are  then  “filtered”  in  a  sense  such  that  dipoles,  a  positive  and  negative 
branch  point  lying  in  adjacent  pixels,  are  cancelled,  as  are  branch  points  which  lie  just 
inside  the  aperture  boundary.  Elsing  the  filtered  branch  point  locations,  the  algorithm 
creates  branch  cuts  to  balance  the  negative  and  positive  residuals.  The  algorithm 
attempts  to  minimize  the  branch  cut  length.  Einally,  it  recursively  unwraps  phase  around 
these  cuts,  yielding  a  path-dependent  phase  retrieval.  Goldstein’s  algorithm  can  fail  if 
branch  cuts  completely  isolate  a  region  of  the  incident  field,  but  this  does  not  seem  to 
occur  in  optical  propagations  [2,4].  The  practical  drawback  to  implementation  of 
Goldstein’s  algorithm  is  that  it  requires  a  wrapped  phase  input,  which  is  generally  not 
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available  in  optical  systems.  Still,  Goldstein’s  performance  provides  a  good  upper  bound 
on  branch  point  tolerant  reconstructor  performance. 

This  chapter  has  presented  the  general  considerations  for  implementing  phase 
reconstructors  in  simulation  and  modeling  atmospheric  turbulence.  Modeling  specifics 
must  also  be  discussed  to  ensure  that  simulation  results  are  reasonable.  These  specifics 
are  given  in  Chapter  4. 
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IV.  Simulation 


Wave-optics  simulations  generally  model  turbulence  via  thin  phase  screens  and 
perform  propagations  using  scalar  diffraction  theory.  “Thin”  phase  screens  are  those 
with  thickness  much  smaller  than  the  propagation  distance  following  the  screen  [29]. 

The  transmittance  of  each  screen  is  given  as 

=  (38) 


where  ¥{x,y)  is  a  zero-mean,  Gaussian  random  phase  with  a  correlation  function  that 
matches  those  of  atmospheric  phase  perturbations.  The  effects  of  this  phase  are  applied 
to  incident  fields  via  Equation  (21). 

In  simulation,  field  propagations  over  long  distances  are  conducted  using  the 
Fresnel  diffraction  integral. 


U{x,y)  =  —e^^^  ^ 


j  \\u{^,V)e 


d^drj 


which  is  clearly  a  scaled  2-D  Fourier  transform  of  the  field  times  a  quadratic  phase  term 
[30].  For  short  distances,  the  angular  spectrum  propagator 


f/(x,y,z)=  j  exp 


xcivcy^ja  -i-|3  jexp  j2n 


«... 


is  used,  where  A 


is  the  angular  spectrum  of  U{x,y,z)  at  z  =  0  and  the 


parameters  a  and  |3  represent  direction  cosines  which  measure  the  angle  of  propagation 


from  the  x  and  y  axes  respectively  [30].  Since  both  propagators  are  based  on  the  Fourier 
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transform,  grid  sampling  and  spacing  is  important  to  proper  modeling.  The  proper  ehoice 
of  these  parameters  depends  on  optieal  geometry,  turbulence,  and  the  partieular 
simulation  method. 


4.1  General  Simulation  Details 


The  simulations  eonsist  of  three  portions,  atmospherie  propagation  of  the  beacon 
field,  phase  reeonstruetion,  and  back  propagation  of  a  corrected  uniform  field  to  the 
beacon  plane.  Propagations  are  conducted  using  WaveTrain™  [31],  an  advanced  wave 
optics  tool  published  by  MZA  Assoeiates  Corp.  The  propagation  geometry  is  shown  in 
Figure  4.  The  beaeon  is  modeled  either  as  an  ideal  point  souree  or  as  an  extended 
beacon.  The  beacon  field  propagates  through  a  series  of  10  evenly  spaced  phase  sereens 
to  a  transmit/receive  aperture.  The  choiee  of  phase  sereens  was  based  on  a 
reeommendation  from  the  sponsor.  Phase  sereens  are  generated  to  give  the  desired 


overall  turbulenee  strength  along  the  path. 
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Figure  4;  Propagation  geometry  for  wave-opties  simulations 


In  order  to  avoid  aliasing  of  the  beacon,  energy  that  does  not  propagate  to  the 
neighborhood  of  the  aperture  is  attenuated.  Simulations  are  condueted  to  model  l-5ms  of 
real-time  behavior  with  field  data  recorded  every  1ms  of  real  time.  To  establish  a 
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statistical  sample,  20  independent  realizations  are  eondueted  for  eaeh  eombination  of 
eontrolled  parameters.  The  ineident  field  at  the  reeeiver  is  measured  after  eaeh 
propagation  using  a  simple  field  sensor  with  resolution  2^x2^  pixels.  Following  eaeh 
beaeon  field  propagation,  phase  is  reeonstrueted  from  the  reeeived  field  using  the  various 
reeonstruetors. 

Baek  propagations  are  eondueted  by  eonverting  the  reeonstrueted  phase  outputs  to 
optieal  path  delays  in  meters.  This  is  done  by  the  following  equation: 

OPD  =  ^xx^^{x,y)  (41) 

where  t|/^  {x,y)  is  the  grid  of  reeonstrueted  phase. 

These  delays  are  then  applied  to  a  uniform  plane  wave,  simulating  the  effeets  of  a 
deformable  mirror,  and  the  resulting  wave  is  back-propagated  from  the  transmit/reeeive 
aperture  to  the  beaeon  plane.  The  on-axis  intensity  is  then  measured  and  the  point  Strehl 
ealculated  as  follows: 

^  ,,  Ineident  On- Axis  Intensity 

Strehl  = - - - .  (42) 

Diffraetion  Limited  On-Axis  Intensity 

4.2  Reeonstruetor  Algorithms 

Three  reeonstruetor  algorithms  are  used  in  the  simulations.  The  hidden  phase 
reeonstruetor  and  Goldstein’s  algorithm  are  branch  point  tolerant.  Least  square  is  also 
implemented  to  provide  a  performanee  eomparison  to  eurrent  systems. 

4.2.1  Least  Squares  Reeonstruetor 

The  Least  Squares  reeonstruetor  diseussed  in  Chapter  2  is  implemented  using  two 
different  eodes.  When  wrapped  phase  input  is  available,  a  diserete  eosine  transform 
algorithm  written  by  Ghiglia  and  Pritt  [7]  is  used.  For  this  eode,  the  wrapped  phase  input 
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must  be  defined  on  a  grid  size  of  (2“+l)x  (2"+l).  Sinee  the  field  sensor  used  is  always  of 
order  2"x2",  a  row  and  eolumn  of  zeros  is  added  to  the  measured  wrapped  phase  before  it 
is  passed  to  the  reeonstruetor.  A  pupil  funetion  is  used  to  limit  the  phase  reeonstruetion 
to  the  region  of  the  aperture.  As  sueh,  the  addition  of  the  zero  row  and  eolumn  has  no 
effeet  on  the  phase  reeonstruetion.  Interfacing  between  WaveTrain  and  the  C  code  is 
done  through  a  Matlab  MEX  system. 

When  wrapped  phase  input  is  not  available,  a  sparse  matrix  based  least  squares 
reconstructor  is  used.  The  reconstructor  calculates  the  matrix  P  given  in  Equation  (1)  and 
applies  sparse  matrix  techniques  to  solve  Equation  (2),  yielding  the  least  squares  solution. 
This  algorithm  requires  only  the  wrapped  phase  differences  derived  from  Equations  (24) 
and  (25)  as  inputs.  This  code  is  accomplished  entirely  in  Matlab  and  was  written  by  Jeff 
Barchers  and  Brent  Ellerbrook  from  the  Air  Eorce  Research  Eab  Starfire  Optical  Range 
[32].  Sample  reconstructions  using  both  algorithms  on  the  same  data  showed  that  the 
methods  were  equivalent. 

4.2.2  Hidden  Phase  Reconstructor 

The  hidden  phase  reconstructor  is  implemented  with  Equation  (36)  in  Matlab. 
Branch  point  locations  are  found  using  Equation  (23),  but  a  cutoff  value  for  branch  point 
detection  is  declared  rather  than  searching  for  exact  values  27i  or  -2%.  Eor  noise-free 
measurements,  a  cutoff  of  ±0.01  provides  for  full  branch  point  detection  with  no  false 
detections.  The  hidden  phase  result  calculated  is  then  summed  with  the  output  of  the 
least  squares  reconstructor  in  use  and  converted  to  optical  path  delay  in  WaveTrain. 

While  simple  to  implement,  this  code  requires  looping  over  detected  branch  points  and 
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can  be  quite  slow  in  strong  turbulence.  While  this  is  somewhat  unimportant  for  a 
simulation  study,  some  optimization  would  be  desirable  in  a  practical  implementation. 

4.2.3  Goldstein’s  Algorithm 

Goldstein’s  algorithm  is  implemented  using  C  code  from  Ghiglia  and  Pritt  [7]. 
This  code  aecepts  wrapped  phase  as  an  input.  A  pupil  aperture  mask  is  also  specified  to 
prevent  reeonstruction  outside  of  the  aperture.  Goldstein’s  algorithm  calculates  the  phase 
differenees  for  branch  point  detection  similar  to  that  used  in  the  hidden  phase  proeess 
using  a  eutoff  value.  The  default  cutoff  value  was  ±0.01,  but  the  eode  has  been  modified 
so  that  the  cutoff  value  can  be  user  speeified  at  run  time. 

Once  the  branch  point  locations  are  determined,  the  code  searches  for  any  positive 
and  negative  branch  points  on  adjacent  pixels.  This  is  ealled  dipole  removal.  These 
pixels  are  removed  from  consideration  and  eonsidered  balanced,  as  are  pixels  which  lie 
adjaeent  to  the  pupil  mask.  Braneh  cuts  are  then  placed  using  the  method  deseribed  in 
Chapter  2.  The  phase  is  then  unwrapped  around  the  cuts  and  saved  as  an  output. 

When  wrapped  phase  is  not  available  as  an  input,  Goldstein’s  algorithm  eannot  be 
used.  The  original  Matlab  interfaee  to  the  Goldstein  C  code  was  first  used  by 
Roggemann  [4]  and  has  been  substantially  modified  and  expanded  for  this  study. 

4.3  Propagation  Geometries 

Two  propagation  geometries  are  used  during  this  study.  The  first  is  patterned 
after  Roggemann’ s  Goldstein  study  [4]  to  allow  for  the  calibration  of  our  simulation. 

This  baseline  provides  a  match  to  peer-reviewed,  published  data,  helping  to  ensure  the 
basic  simulation  issues  are  handled  correctly.  The  second  geometry  is  based  on  the  Air 
Force  Research  Laboratory  North  Oscura  Peak  (NOP)  optical  range  in  New  Mexico  [33]. 
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4.3.1  Baseline  Geometry 

The  first  geometry  is  based  on  Roggemann’s  Goldstein  study  [4]  to  allow  for  the 
ealibration  of  our  simulation.  A  0.5m  telescope  is  used  over  a  100km  propagation  path. 
The  atmospheric  grid  spacing  is  set  to  1cm  with  a  point  source  beacon  of  wavelength 
0.987  pm.  The  field  sensor  and  reconstructor  dimensions  are  256x256  for  most  of  the 
baseline  geometry  studies. 

The  validity  of  our  atmospheric  phase  screens  is  tested  by  calculating  the  log- 
amplitude  variance  in  the  receiver  aperture  against  the  theoretical  Rytov  variance  for  each 
turbulence  strength.  It  is  well  documented  that  the  log-amplitude  variance  saturates  in 
real-world  conditions  and  simulation  [28].  As  Figure  5  shows,  the  experimental  data 
provides  a  good  fit  to  theory  up  to  the  saturation  region.  Figure  5  is  also  a  close  match  to 
Figure  4  in  Roggemann’s  study  [4].  Since  Roggemann  used  a  completely  different  wave- 
optics  simulation  for  his  propagations,  this  similarity  shows  that  WaveTrain  should 
provide  nearly  equivalent  results. 

Three  studies  were  conducted  using  the  baseline  geometry  -  noise-free/zero  time- 
delay,  noise-free/fmite  time  delay,  and  noise-corrupted/zero  time  delay.  Each  study  used 
a  point  source  beacon.  The  case  of  noise-corrupted/fmite  time  delay  was  not  considered 
since  the  two  effects  are  independent.  System  performance  under  a  combination  of  noise 
and  time  delay  can  be  reasonably  inferred  from  the  individual  results  of  the  noise- 
corrupted  and  time  delayed  studies.  The  noise-free,  zero  time-delay  study  was  used  to 
establish  a  performance  baseline  for  each  reconstructor  and  compare  our  simulation 
performance  to  Roggemann’s  results.  It  also  represents  the  ideal  AO  system. 
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The  noise-free,  finite  time-delay  shows  the  impact  of  finite  time-delayed 
correction  on  system  performance.  For  this  study,  the  phase  reconstruction  is  calculated 
using  the  initial  beacon  propagation  and  then  held  for  4ms  on  the  back  propagation.  A 
constant  wind  was  used  to  change  the  turbulence  in  the  propagation  path.  The 
atmospheric  shift  is  easily  calculated  by  multiplying  the  wind  velocity  by  the  correction 
delay  time.  Four  wind  speeds  are  used;  25,  50,  100,  and  150  knots.  These  were  chosen 
to  provide  a  good  sample  of  shift  distances.  The  possible  shifts  are  1-30  centimeters  for 
time  delays  of  1  -  4  ms. 

Finally,  a  noise  study  was  conducted  using  the  hidden  phase  processor.  For  this 
study,  white  Gaussian  noise  of  various  strengths  is  added  to  the  phase  difference 
measurements  prior  to  branch  point  detection  and  least  squares/hidden  phase 
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reconstruction.  To  account  for  this  noise,  various  cutoff  values  are  used  to  declare  branch 
points  in  Equation  (23)  in  order  to  see  their  impact  on  performance. 

4.3.2  North  Oscura  Peak  (NOP)  Geometry 

For  this  geometry,  the  receiver  is  modeled  as  a  0.75  m  telescope  focused  on  the 
beacon  source.  An  ideal  field  sensor  measures  the  received  field  for  phase 
reconstruction.  A  beacon  wavelength  of  1 .0  pm  is  used.  The  total  propagation  distance 
is  50  km  with  an  atmospheric  grid  sample  spacing  of  2  cm. 

Three  studies  are  conducted  for  the  NOP  geometry  -  ideal  point  source, 
incoherent  extended  beacon,  and  coherent  extended  beacon.  Each  of  these  studies  uses 
noise-free,  zero  time-delay  corrections  to  emphasize  the  beacon  contribution  to  overall 
performance.  The  ideal  point  source  study  is  used  to  form  the  ideal  baseline  for  NOP, 
similar  to  that  done  for  the  first  geometry.  Details  of  the  extended  beacon  studies  will 
follow  in  Section  4.6. 

4.4  Choice  of  Modeling  Parameters 

The  choice  of  correct  modeling  parameters  is  critical  in  wave-optic  simulations. 
Since  atmospheric  propagations  are  conducted  through  transforms,  attention  must  be 
given  to  spatial  sampling  to  avoid  aliasing.  The  field  sensor  model  is  also  important,  as  it 
tends  to  drive  the  computational  requirements  for  reconstruction. 

4.4.1  Ideal  Field  Sensor  Resolution 

The  computational  time  and  disk  storage  requirements  for  phase  reconstruction 
are  dependent  primarily  on  the  pixel  count  of  the  incident  wrapped  phase.  As  such,  it  is 
desirable  to  use  the  smallest  grid  resolution  possible,  which  does  not  largely  impact 
overall  performance.  Goldstein’s  algorithm  requires  that  grid  resolutions  be  powers  of  2 


34 


while  one  of  the  least  squares  algorithms  implemented  uses  the  discrete  cosine  transform 
and  requires  that  the  resolution  be  an  even  power  of  2  plus  1  pixel  (2°+l). 

Following  the  baseline  studies  listed  in  Section  4.4.1,  an  additional  study  was 
conducted  using  the  baseline  geometry  to  determine  Goldstein  Strehl  performance  vs. 
field  sensor  resolution.  The  choice  of  field  sensor  resolutions  was  ad-hoc,  but 
encompassed  that  of  current  wave-front  sensors  and  some  extreme  values.  The  six 
implemented  NxN  resolutions  were  N  =  [8,  16,  32,  64,  128,  256]  which  are  powers  of 
two  and  can  be  directly  used  with  Goldstein’s  algorithm.  Figure  6  shows  the  relative 
performance  for  each  N  using  a  point  source  beacon.  As  can  be  clearly  seen  in  the  figure, 
there  is  only  a  minor  difference  between  the  N=256  and  N=128  pixel  resolutions, 
especially  in  the  Rytov  range  below  6  where  the  NOP  simulations  lie.  Based  on  this 
study,  it  was  decided  to  use  a  128x128  field  sensor  for  Goldstein’s  algorithm  with  the 
NOP  geometry.  It  was  also  decided  to  use  this  resolution  for  least  squares  reconstruction 
to  maintain  consistency  in  simulation  geometry  and  data  handling.  Least  squares 
reconstruction  is  much  less  sensitive  to  sensor  resolution. 

4.4.2  Choice  of  Atmospheric  Propagation  Grid  Resolution 

The  extent  of  the  atmospheric  propagation  grid  must  be  sufficient  to  avoid  wrap¬ 
around  during  the  FFT  propagations.  The  individual  pixel  widths  must  also  be  small 
enough  to  capture  sufficient  detail  considering  the  diffraction  limit  of  our  system. 
WaveTrain’s  resolution  is  specified  by  the  pixel  width  Ax  and  pixel  grid  size  NxN.  The 
system  diffraction  limit  is  given  as 

Xz 

D.L.  =  —  (43) 

D 
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Figure  6;  Strehl  vs.  Rytov  Variance  for  Various  Field  Sensor  Resolutions 


where  X  is  the  wavelength,  z  the  propagation  distance,  and  D  the  diameter  of  our  optic. 
Using  >.  =  1  pm,  z  =  50  km,  and  D  =  0.75  m,  which  are  the  general  parameters  for  the 
NOP  geometry,  D.L.  is  0.067  m  or  roughly  7  cm.  This  is  the  largest  pixel  width  we 
should  consider  for  NOP  simulations. 

We  also  need  to  consider  turbulence  sampling,  which  is  related  to  the  Fried 
parameter  tq  .  For  a  constant  C„  along  the  propagation  path,  this  is  expressed  as 

r,=  2.1(1  (44) 

where  z,  is  the  distance  between  phase  screens.  The  smallest  ro  in  the  NOP  simulations  is 

about  7.5  cm.  Based  on  recommendations  from  the  thesis  sponsor,  2  cm  was  used  for  Ax. 
This  choice  ensures  that  there  are  multiple  samples  across  the  smallest  expected  ro.  N 
must  then  be  chosen  large  enough  to  prevent  wrap-around.  The  WaveTrain 
documentation  [31]  suggests  Nmin  be  chosen  using  the  following  equation: 
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(45) 


N„ 


iXz 


which  gives  an  Nmin  of  250  for  NOP.  Since  N  must  be  a  power  of  2,  256  is  the  smallest 
possible  value. 

Figure  7  shows  the  baseline  runs  for  Goldstein’s  algorithm  using  NOP  geometry 
with  N  =  256  and  Ax  =  2  cm.  It  is  clear  in  the  plot  that  some  type  of  modeling  issue 
occurs  for  Cl  larger  than  5x10'^^  This  is  evidenced  by  the  upturns  in  the  Strehl 
performance  where  none  are  expected.  This  error  is  due  to  insufficient  atmospheric  grid 
size  for  these  turbulence  strengths  and  wraparound  has  occurred  during  propagation. 
Based  on  this  information,  atmospheric  grid  size  was  increased  to  N  =  512  beginning  at 

=  5x  10'  m'  .  A  constant  C„  was  used  for  simplicity,  but  this  research  could  be 

2 

easily  repeated  for  specific  C„  models  of  interest. 


Figure  7:  Strehl  vs.  Field  Sensor  Resolution  for  NOP  Geometry 
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4,5  Extended  Beacon  Modeling  Issues 

Wave-optics  simulations  are  performed  based  on  the  theory  of  coherent  field 
propagation.  Reconstructor  theory,  on  the  other  hand,  is  based  on  an  assumption  of  an 
incoherent  beacon  [21],  The  point  source  beacon  provides  a  special  case  because  it  is  a 
coherent  source,  but  by  lacking  finite  extent,  it  appears  incoherent  to  an  optical  system. 
Although  true  point  sources  do  not  exist  in  nature,  beacons  with  finite  extent  less  than  the 
diffraction  limit  of  an  optical  system  appear  as  point  sources,  hence  the  success  in  the  use 
of  adaptive  optics  in  astronomical  imaging. 

In  the  case  of  Airborne  Laser  and  other  systems  which  must  use  artificial  beacons, 
the  beacon  extent  may  exceed  the  system  diffraction  limit.  Such  extended  beacons  are 
generally  considered  incoherent  due  to  various  factors.  Most  physical  surfaces  contain 
microscale  roughness  on  the  order  of  optical  wavelengths,  thus  causing  random  spatial 
coherent  speckle  on  reflection.  Also,  large-scale  target  depth,  which  exceeds  the 
coherence  length  of  the  illuminator  laser,  reduces  the  coherence  of  the  reflected  beacon. 
WaveTrain  does  not  currently  support  such  a  detailed  propagation  and  reflectance  model 
and  other  methods  must  be  used  to  approximate  an  incoherent  extended  beacon. 

4.5.1  Speckle  Method  for  an  Incoherent  Extended  Beacon 

Modeling  extended  incoherent  beacons  using  wave-optics  simulation  is 
particularly  difficult.  The  simulated  beacon  source  must  be  coherent  in  order  to  perform 
propagations.  Incoherence  must  be  modeled  at  the  receiver  using  the  results  of  coherent 
beacon  propagations. 
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For  this  research,  two  methods  were  examined  for  simulating  an  incoherent 
extended  beacon.  The  use  of  multiple  speckle  fields  was  the  first  technique.  This 
technique  was  implemented  in  WaveTrain  using  an  incoherent  reflector  model.  The 
incoherent  reflector  creates  speckle  realizations  using  the  following  three  steps.  First,  the 
user  specifies  an  amplitude  reflectance  map.  WaveTrain  then  reflects  a  user-defined  field 
from  the  intensity  map  to  create  a  new  optical  field.  For  this  thesis,  the  user-defined 
illuminator  field  is  always  a  uniform  coherent  source.  Prior  to  any  propagation,  the  phase 
of  the  new  field  is  replaced  with  a  zero-mean  random  phase  which  is  uniform  on  -n  to  n. 
A  spatial  fdter  is  applied  as  a  final  step  before  propagation  to  avoid  wraparound.  The 
field  is  then  propagated  to  the  receiver  as  a  single  speckle  realization.  The  process  is 
repeated  for  each  speckle  realization. 

Reflectance  maps  are  created  in  Matlab  using  a  Gaussian  profile.  Care  is  taken  to 
ensure  that  the  reflectance  is  near  zero  at  the  map  edges  to  prevent  ringing  or  aliasing. 

For  this  study,  beacon  extent  always  refers  to  the  2xa  value  for  a  standard  Gaussian 
profde 


F  = 


V^(7 


-exp 


1  f 


(46) 


with  r  equal  to  the  radius  from  the  x-y  origin  of  the  beacon.  The  beacon  is  always 
centered  on  the  optical  axis. 

The  number  of  speckle  realizations  required  should  depend  on  the  beacon  size  and 
turbulence  strength.  We  initially  examined  beacon  extents  out  to  0.5m.  Since  coherent 
speckle  decreases  as  turbulence  strength  increases  due  to  the  mixing  in  the  propagation 
path,  vacuum  propagation  represents  the  worst  case  scenario. 
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To  determine  the  requirement  on  speekle  realizations,  a  study  was  eondueted 
using  vaeuum  propagation  for  four  beaeon  sizes.  First,  250  independent  speekle 
realizations  were  generated  for  eaeh  beaeon  size.  The  ineident  phase  was  ealeulated 
using  Equation  (20).  This  result  is  always  wrapped  on  -n  to  n  due  to  the  tangent  inverse 
operation. 

Averages  of  the  wrapped  phase  at  the  field  sensor  were  then  taken  over  5  to  250 
realizations  in  steps  of  5.  Goldstein’s  algorithm  was  then  used  to  form  a  reeonstrueted 
phase  for  eaeh  wrapped  phase  average.  These  reeonstrueted  phases  were  then  used  as 
optieal  path  delays  in  WaveTrain  and  baek-propagated  to  the  beaeon  plane.  Figure  8 
shows  the  resulting  Strehl  ratios  as  a  funetion  of  phase  realizations  averaged  and  beaeon 
size.  We  see  that  as  the  beaeon  size  inereases,  more  independent  phase  realizations  are 
needed  for  the  reeonstrueted  average  phase  to  aehieve  95%  of  the  uneorreeted 
performanee  (Strehl  of  1.0). 

Based  on  the  results  in  Figure  8,  twenty  independent  point  Strehl  averages  were 
obtained  for  speekle  averages  of  10  to  100  in  steps  of  10.  The  goal  for  this  study  was  to 
find  the  number  of  speekle  averages  required  to  aehieve  95%  of  the  ideal  Strehl 
performanee  plus  1  standard  deviation.  These  results  indieated  that  60  speekles  appeared 
to  be  appropriate  for  modeling  the  ineoherent  extended  beaeon. 

Unfortunately,  this  teehnique  fails  to  provide  meaningful  results  when  applied  to 
turbulenee  propagations.  An  examination  of  the  averaged  ineident  phases  reveals  that  the 
mean  of  the  average  phase  approaehes  zero  as  the  number  of  phase  realizations  averaged 
inereases.  This  is  due  to  the  ineoherent  refleetor  using  a  zero-mean  phase  for  eaeh 
speekle  realization.  In  vaeuum,  the  eorreet  phase  reeonstruetion  should  provide  a  zero 


40 


path-delay,  and  hence  the  technique  showed  consistent  improvements  in  Strehl  ratio  as 
the  number  of  phase  realizations  averaged  increased. 


Figure  8;  Strehl  vs.  Number  of  Speckle  Averages  in  Vacuum  (Method  1) 


In  turbulence,  the  receiver  incident  phase  is  the  wrapped  sum  of  a  non-zero  mean 
atmospheric  phase  (a  fixed  random  variable  for  each  speckle  average)  and  the  zero-mean 
speckle  phases.  The  wrapping  function  on  this  sum,  which  is  a  non-zero  mean  Gaussian 
random  variable,  maps  the  wrapped  phase  average  to  zero  mean.  Thus,  as  the  number  of 
speckle  averages  increases,  the  average  phase  converges  toward  zero  and  all  information 
about  the  atmospheric  turbulence  is  averaged.  The  resulting  phase  reconstruction 
provides  a  near-zero  path-delay  for  all  turbulence  strengths.  This  results  in  performance 
within  2%  of  the  uncompensated  performance.  Clearly,  another  method  must  be  used  to 
model  the  beacon. 

The  proper  technique  for  modeling  the  incoherent  extended  beacon  is  to  average 
the  wrapped  phase  differences.  A""  and  A^' ,  given  by  Equations  (24)-(25).  This  technique 
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is  similar  to  the  integration  whieh  oeeurs  in  a  praetieal  wave  front  sensor  over  a 
measurement  sample  time.  Unfortunately,  Goldstein’s  algorithm  eannot  be  used  with  this 
teehnique,  as  the  wrapped  phase  input  is  unavailable.  Deriving  wrapped  phase  from 
phase  differenees  is  a  phase  retrieval  proeess  not  unlike  finding  the  unwrapped  phase. 
Least  squares  ean  be  used  to  estimate  a  wrapped  phase,  but  by  nature,  it  eliminates 
braneh  point  information  eritieal  to  Goldstein’s  reeonstruetor.  As  sueh,  least  squares  and 
hidden  phase  are  used  for  the  incoherent  beacon  study  since  they  are  formulated  for  phase 
differences. 

The  choice  of  speckle  realization  was  once  again  determined  by  conducting 
vacuum  propagations.  The  desired  performance  threshold  was  0.95  Strehl  and  no  false 
branch  point  detections  using  a  0.2  cutoff  value  for  declaring  if  a  phase  residual  indicated 
a  branch  point.  Three  beacon  sizes  were  considered  -  4cm,  7.5cm,  and  16cm.  These  are 
essentially  0.5,  1,  and  2x  the  diffraction  limit  of  NOP.  The  vacuum  study  determined 
that  the  4cm  beacon  required  10  speckle  averages,  the  7.5cm,  15,  and  the  16cm,  30  to 
achieve  the  performance  threshold.  As  the  results  will  show,  this  technique  was  also 
appropriate  for  turbulence  propagations. 

4.5.2  Multiple  Point  Source  Model  for  a  Coherent  Extended  Beacon 

Multiple  point  sources  were  used  in  the  second  extended  beacon  model.  The 
desire  was  to  obtain  a  model  useable  with  Goldstein’s  algorithm.  For  this  method,  an 
array  of  independent  point  sources  was  created  having  a  desired  spatial  extent  and 
density.  The  intensity  of  each  point  source  was  set  according  to  its  grid  location  to  create 
a  Gaussian  intensity  profile.  The  point  source  fields  are  individually  propagated  to  the 
receiver  and  then  summed  together  coherently.  This  forms  a  coherent  Gaussian  beacon 
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model.  Consideration  was  given  to  using  an  existing  eoherent  Gaussian  source  in 
WaveTrain,  but  initial  attempts  to  set  up  this  source  gave  inconsistent  results. 

A  study  was  conducted  to  determine  the  number  and  density  of  point  sources 
required  to  achieve  consistent  performance.  For  the  study,  the  two  parameters  of  interest 
were  the  spacing  between  point  sources  and  the  cut-off  value  in  meters  for  the  total  non¬ 
zero  beacon  extent.  The  study  was  conducted  for  vacuum  propagations,  which  represent 
the  worst-case  scenario  since  the  system  resolution  is  at  its  highest. 

It  was  found  that  both  least  squares  and  Goldstein’s  algorithms  were  very 
sensitive  to  the  “smoothness”  of  the  Gaussian  cut-off  and  the  “quality”  of  the  Gaussian 
profile.  By  experiment,  it  was  found  that  a  cut-off  of  3a  combined  with  a  point  source 
every  3  cm  in  the  grid  gave  best  results.  The  Gaussian  intensity  at  3a  is  only  0.4%  of  the 
starting  value  -  thus  3  a  gives  a  very  smooth  cut-off  The  diffraction  limited  resolution 
for  the  NOP  geometry  is  ^1  cm,  so  each  point  source  is  roughly  half  a  diffraction  limited 
spot  size  from  its  nearest  neighbor.  Thus,  the  point  sources  are  not  individually 
resolvable  by  the  system,  which  helps  provide  a  good  Gaussian  approximation. 

After  propagating  this  type  of  beacon  in  vacuum,  we  hope  to  see  a  nice  Gaussian 
intensity  and  a  uniform  phase  similar  to  that  in  Figure  1 .  The  corresponding  correction 
should  be  zero  and  we  should  get  a  Strehl  of  '-1 .0  after  back  propagation.  This  is  not 
what  occurs. 

Figure  9  shows  the  performance  for  five  beacon  sizes  using  vacuum  propagation. 
Clearly,  system  performance  starts  to  degrade  quickly  for  beacons  larger  than  the 
diffraction  limit.  This  is  due  to  coherent  speckle  caused  by  the  coherent  beacon.  The 
coherent  speckle  induces  amplitude  perturbations  in  the  field,  resulting  in  branch  points. 
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The  wave-front  sensor  measures  these  phase  effeets  as  it  would  for  those  eaused  by 
turbulenee.  The  system  has  no  way  to  distinguish  between  field  perturbations  eaused  by 
turbulenee  and  field  perturbations  eaused  by  eoherent  speekle.  Sinee  we  are  only 
interested  in  measuring  the  turbulenee  effeets,  the  eoherent  speekle  ereates  a  form  of 
measurement  error  that  grows  with  beaeon  size.  This  error  ereates  errors  in  the  phase 
reeonstruetion  and  reduees  Strehl  performanee.  It  is  elear  from  Figure  9  that  studying 
beacon  sizes  above  0.5  m  holds  little  interest,  as  the  vacuum  performance  is  terrible. 


Figure  9:  Strehl  vs.  Beacon  Diameter  in  Vacuum  for  5  Extended  Beacons 


Propagating  the  point  source  grids  requires  considerable  computation  time.  Using 
Figure  7  as  justification,  it  was  decided  to  concentrate  the  coherent  extended  beacon 
Study  to  C„  between  UlO"  and  9x10'  m'  to  reduce  the  sample  set.  In  weaker 
turbulence,  the  performance  of  Goldstein  and  least  squares  reconstruction  is  essentially 
the  same  since  branch  points  are  not  generally  present.  The  Rytov  number  at 


44 


=  1 X 10'*^  using  the  NOP  geometry  is  0.0436,  slightly  below  the  value  of  the 
expected  onset  of  branch  points  above  a  Rytov  of  0.05. 

By  establishing  proper  modeling  parameters,  we  can  have  confidence  in  our 
simulation  results.  The  extended  beacon  models  will  provide  new  insights  into  branch 
point  reconstructor  performance.  The  choice  of  propagation  geometries  allows  us  to 
verify  our  ideal  simulations  against  published  work,  while  allowing  us  to  scale  our  non¬ 
ideal  performance  to  systems  of  interest  like  the  Airborne  Laser. 
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V.  Results 


This  chapter  presents  results  from  eaeh  of  the  three  major  studies.  The  baseline 
studies  are  useful  to  establish  upper  bounds  on  ideal  performanee.  They  also  eonsider  the 
effeets  of  time-delays  in  applying  optieal  eorreetions  and  give  some  insight  into  the 
problem  of  finding  braneh  points  with  noise-eorrupted  phase  measurements.  The 
extended  beacon  studies  establish  performanee  bounds  for  two  important  limiting  oases, 
full  beaoon  ooherenoe  and  full  beaoon  inooherenoe. 

5.1  Baseline  Geometry  Results 

The  baseline  geometry  is  important  to  establish  the  validity  of  our  simulations  and 
provide  an  ideal  upper  bound  on  performanee.  The  ideal  performanee  for  zero-time- 
delay  eorreetions  and  noise-free  measurements  was  very  olose  to  published  results  [4]. 
Non-ideal  soenarios  with  time  delay  and  measurement  noise  give  insight  into  the 
robustness  of  the  various  reoonstruotors. 

5.1.1  Ideal  Performance 

Performanee  for  noise-free,  zero-time-delay  eorreetions  is  shown  in  Figure  10  and 
Figure  1 1 .  Goldstein’s  algorithm  provides  the  best  performanee  in  strong  turbulenee 
(Rytov  >  0.5)  with  a  nearly  2  orders  of  magnitude  improvement  over  uneompensated 
performanee.  The  hidden  phase  algorithm  has  elear  problems  at  Rytov  <1.3,  but  does 
somewhat  better  than  least  squares  above  a  Rytov  of  1 .3.  The  least  squares  algorithm 
matehes  Goldstein  performanee  out  to  about  0.5  Rytov,  beyond  whieh  braneh  points 
dominate  reeonstruetor  performanee. 
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Strehl 


Figure  1 1 :  Strehl  vs.  Rytov  for  Noise-free  Zero  Time  Delay  Baseline  Geometry 
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5.1.2  Time-delayed  Correction  Performance 

Figure  12  through  Figure  14  illustrate  the  performanee  of  eaeh  reeonstruetor  for 
finite  eorreetion  time  delays  with  the  25  knot  wind.  For  eomparison,  the  performanee 
versus  time  delay  for  the  150  knot  wind  is  shown  in  Figure  15-  Figure  17.  Sinee  a 
eonstant  wind  was  used  to  move  the  atmosphere,  we  would  expeet  that  the  performanee 
degradation  should  seale  linearly  with  wind  speed.  Based  on  this  model,  we  ean  also 
approximate  the  eorrelation  distanee  of  the  atmosphere  by  multiplying  the  wind  veloeity 
by  the  time  delay  required  for  compensated  performance  to  approximate  uncompensated. 
Experimental  correlation  times  and  distances  for  our  simulations  are  listed  in  Tables  1 
and  2.  Times  were  rounded  to  the  nearest  1ms  due  to  our  time  discretization. 

An  examination  of  the  results  shows  that  performance  is  tied  to  shift  distance  as 
expected.  The  shift  distance  for  decorrelation  was  similar  for  all  three  reconstructors 
with  the  branch  point  tolerant  reconstructors  slightly  more  sensitive  to  shift.  This  is  most 
likely  due  to  the  localized  nature  of  branch  points.  There  is  also  a  clear  difference  in  the 
correlation  distance  for  cl  >  7x10'^^  For  >  7x10'^^  the  correlation 
distance  was -O.lrn  on  average.  Correlation  distance  for  C„  <7x10'  m'  was  about 
0.2m,  except  for  the  hidden  phase  algorithm  which  decorrelated  somewhat  faster  than  the 
others.  We  note  these  values  are  on  the  order  of  the  Fried  parameter,  ro,  for  these 
turbulence  regions.  For  long  time  delays  and  strong  turbulence,  corrected  performance 
approaches  uncompensated  performance,  indicating  that  the  corrections  are  decorrelated 
from  the  atmosphere.  The  data  also  indicates  that  applying  decorrelated  corrections  can 
result  in  Strehl  ratios  worse  than  uncompensated,  but  rarely  much  worse. 
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Figure  12:  Strehl  vs.  Cl  and  Time  Delay  for  Goldstein  with  a  25  Knot  Wind 
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Figure  14:  Strehl  vs.  C„  and  Time  Delay  for  Least  Squares  with  a  25  Knot  Wind 


Figure  15:  Strehl  vs.  Cl  and  Time  Delay  for  Goldstein  with  a  150  Knot  Wind 
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Figure  16;  Strehl  vs.  C„  and  Time  Delay  for  Hidden  Phase  with  a  150  Knot  Wind 


2 

Figure  17:  Strehl  vs.  C„  and  Time  Delay  for  Least  Squares  with  a  150  Knot  Wind 
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Table  1:  Correlation  times/distances  for  >  7x10'^^ 


25Kts 

50Kts 

lOOKts 

150Kts 

Goldstein 

>4ms 

3ms/  0.08m 

2ms/  0.102m 

2ms/  0.152m 

Hidden  Phase 

>4ms 

3ms/  0.08m 

2ms/  0.102m 

2ms/  0.152m 

Least  Squares 

>4ms 

4ms/  0.102m 

3ms/  0.152m 

2ms/  0.152m 

Table  2;  Correlation  times/distances  for  <  7x10'^^ 


25Kts 

50Kts 

lOOKts 

150Kts 

Goldstein 

>4ms 

>4ms 

4ms/  0.203m 

3ms/  0.229m 

Hidden  Phase 

>4ms 

4ms/  0.102m 

3ms/  0.152m 

2ms/  0.229m 

Least  Squares 

>4ms 

>4ms 

4ms/  0.203m 

3ms/  0.229m 

5.1.3  Noise  Study  Results 

The  noise  study  indieates  that  the  hidden  phase  algorithm  can  still  outperform 
simple  least  squares  reconstruction  in  the  presence  of  measurement  noise  if  the  branch 
point  cutoff  is  chosen  correctly.  Figure  18  is  a  plot  of  hidden  phase  performance  for  four 
noise  strengths  and  ten  cutoff  values.  The  cutoff  values  have  been  normalized,  where  a 
cutoff  of  1  represents  a  value  of  27i  in  the  circulation  residual.  The  atmospheric 
turbulence  strength  was  set  at  =  8x10'^^  The  Strehl  values  are  an  average  of  20 
realizations  for  each  combination  of  noise  and  cutoff  value.  The  effects  of  choosing  a 
cutoff  too  close  to  zero  or  one  are  seen  as  the  performance  quickly  degrades.  It  is  also 
apparent  that,  as  the  noise  strength  increases,  the  importance  of  the  specific  cutoff  value 
also  grows  due  to  overlap  in  the  residual  distributions. 

Some  thought  was  given  to  deriving  an  analytic  expression  to  determine  the 
optimum  branch  point  cutoff  Maximum  likelihood  techniques  are  often  appropriate  for 
this  class  of  problem,  but  difficulty  arises  for  the  specific  application  of  branch  point 
finding.  Consider  the  strong-turbulence  histogram  of  noise-free  phase  residuals  from  one 
propagation  given  in  Figure  19.  The  residuals  have  been  normalized  such  that  a 
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Figure  18:  Strehl  vs.  Normalized  Branch  Point  Cutoff  Value  for  four  noise  strengths  at 

C„^of8xl0-^^m-^^^ 
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Figure  19:  Histogram  of  Branch  Point  Residuals 
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positive  branch  points  are  indicated  at  1  and  negative  branch  points  at  -1 .  We  see  that 
the  positive  and  negative  branch  points  occur  in  relatively  equal  numbers,  but  much  less 
frequently  than  non-branch  point  pixels.  The  relative  probabilities  of  occurrence  are 
1.5%  each  for  positive  and  negative  branch  points  and  97%  for  non-branch  point  pixels. 

The  addition  of  white  Gaussian  noise  to  the  phase  differences  can  be  interpreted 
as  a  convolution  of  a  Gaussian  function  with  each  of  the  impulses  in  the  histogram. 

Since  non-branch  point  pixels  contain  most  of  the  probability,  the  center  Gaussian  will 
overlap  the  branch  point  distributions  for  most  noise  strengths  of  interest.  As  this  occurs, 
the  maximum  likelihood  decision  is  to  declare  all  pixels  as  non-branch  point  pixels.  As 
seen  in  Figure  18,  this  is  the  wrong  decision,  as  the  cost  of  missing  branch  points  (high 
cutoffs)  is  high.  We  cannot  set  the  cutoff  too  low  either,  because  there  is  also  a  cost  for 
false  branch  points,  as  seen  in  the  strongest  noise  case. 

5.2  NOP  Geometry  Results 

The  North  Oscura  Peak  geometry  is  important  because  results  obtained  with  this 
geometry  are  scalable  to  Airborne  Laser  and  other  systems.  Branch  point  reconstructors 
are  most  desirable  in  strong  turbulence.  Systems  like  ABL,  which  illuminate  their 
targets,  will  often  have  extended  beacons  in  such  turbulence  so  it  is  critical  to  evaluate 
how  branch  point  reconstructors  perform  with  extended  beacons. 

5.2.1  Incoherent  Extended  Beacon  Performance 

Baseline  performance  for  the  NOP  geometry  is  shown  in  Figure  20  and  Figure  21. 
These  plots  were  obtained  by  running  point  source  simulations  over  the  exact  same 
turbulence  used  with  the  incoherent  extended  beacons.  Reconstructions  were  done  with 
phase  slopes  only;  thus  Goldstein’s  algorithm  was  not  used.  The  small  jump  in 
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2  17  2/^ 

performance  at  C„  =  5x10'  m'  is  due  to  the  shift  in  atmospheric  grid  size  from 
256x256  to  512x512. 

Two  methods  were  used  to  specify  branch  points  for  the  hidden  phase  processor  - 
the  normal  circulation  method  and  the  Goldstein’s  “filtered”  approach.  It  is  evident  that 
the  Goldstein  approach  does  help  in  weak  turbulence,  but  in  strong  turbulence  the 
performance  is  much  worse  than  using  standard  branch  point  detection.  In  weak 
turbulence,  the  circulation  method  sometimes  detects  single  branch  points  along  the 
aperture  edge,  which  are  caused  by  edge  effects  and  not  the  turbulence.  These  edge 
branch  points  are  normally  eliminated  in  the  Goldstein  detection  scheme. 

Goldstein’s  algorithm  also  eliminates  branch  point  dipoles.  In  terms  of  branch 
cuts,  a  dipole  would  represent  a  zero  cut  length  since  we  have  a  discrete  unwrapping. 
Since  it  impossible  to  cross  a  zero-length  cut,  Goldstein  eliminates  branch  point  dipoles. 
The  hidden  phase  process  does  not  rely  on  cuts.  Instead  it  calculates  the  least  squares 
reconstruction  error  for  each  pixel  due  to  the  contributions  of  all  branch  points.  In  this 
case,  the  individual  branch  points  in  a  dipole  remain  significant  and  cannot  be  neglected. 
The  number  of  dipoles  increases  with  turbulence  strength,  and  when  these  are  eliminated, 
hidden  phase  performance  converges  on  simple  least  squares. 

The  hidden  phase  performance  in  weak  turbulence  for  the  NOP  geometry  is 
somewhat  better  than  that  obtained  for  the  baseline  geometry.  This  is  possibly  due  to  the 
use  of  a  128x128  field  sensor  grid,  which  maybe  less  sensitive  to  false  branch  point 
detection  at  the  aperture  edges.  The  break-even  point  between  least  squares  and  hidden 

2  17  2/^ 

phase  occurs  at  essentially  the  same  value  of  C„  --8x10'  m'  or  about  0.5  Rytov. 
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Figure  20:  Strehl  vs.  Cl  for  NOP  Point  Souree  Beacon 


Figure  21:  Strehl  vs.  Rytov  for  NOP  Point  Source  Beacon 
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Figure  22  provides  a  plot  of  Strehl  ratio  vs.  LOo  for  least  squares  reeonstruetion. 
In  this  ease  L  is  the  propagation  distanee,  50  km,  and  6o  is  the  isoplanatie  angle.  Thus, 
L6o  represents  the  isoplanatie  pateh.  In  this  figure,  turbulenee  strength  inereases  from 
right  to  left  and  we  see  that  performanee  degrades  as  the  isoplanatie  pateh  deereases. 


Figure  22:  Strehl  vs.  LO^  for  Least  Squares 


Figure  23  shows  least  squares  performanee  for  the  4,  7.5,  and  16em  ineoherent 
beacons.  As  expected,  the  performance  is  degraded  by  finite  beacon  extent.  It  is  also 
interesting  to  note  that  the  relative  performance  degradation  of  the  4  and  7.5cm  beacons 
gets  larger  with  beacon  size  to  a  Rytov  of  '-1 .6,  and  then  the  performance  delta  narrows 
again  for  larger  Rytov.  As  the  turbulence  strength  increases,  the  spatial  correlation 
decreases,  thus  the  extended  beacons  sample  a  region  larger  than  the  isoplanatie  patch. 
This  causes  errors  in  least  squares  reconstruction  -  with  the  error  increasing  with  the 
proportion  of  beacon  outside  the  isoplanatie  angle.  As  the  16cm  beacon  line  shows,  this 
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error  tends  to  saturate  above  a  eertain  size  as  the  16em  beaeon  performanee  eonverges 
slightly  on  the  8em  performanee  at  Rytov  =  1 .6.  The  eonvergenee  of  eaeh  of  the  beaeons 
in  strong  turbulenee  is  due  to  the  dominanee  of  braneh  point  effeets.  Still,  least  squares 
reeonstruetion  does  provide  a  nearly  5x  improvement  over  the  uneompensated  ease. 

A  major  problem  was  diseovered  with  using  braneh  point  deteetion  algorithms  on 
ineoherent  beaeon  sourees.  This  is  shown  in  Figure  24,  a  plot  of  reeonstruetor 
performanee  for  the  4em  ineoherent  beaeon.  The  beaeon  ineoherenee  has  removed 
braneh  point  information  from  the  eireulation  residuals  through  the  phase  averaging 
proeess.  Without  information  on  braneh  point  loeations,  the  hidden  phase  proeessor 
reverts  to  a  least  squares  reeonstruetor.  The  same  result  oeeurred  for  the  larger  beaeons. 

Examination  of  the  eireulation  values  for  eaeh  extended  beaeon  size  revealed  that 
they  were  on  the  order  of  ±0.2  in  most  oases.  Peak  residuals  did  not  oorrespond  to  true 
braneh  point  loeations  with  any  regularity.  This  indioates  a  fundamental  diffioulty  in 
applying  braneh  point  deteetion  to  ineoherent  extended  beaeons.  As  will  be  shown  in 
Seotion  5.2.2,  braneh  point  information  was  maintained  with  the  ooherent  beaeons. 

Based  on  this  result,  the  loss  of  braneh  point  information  with  the  ineoherent  beaeon  is 
likely  due  to  the  phase  averaging  proeess  whioh  we  simulated  through  speekle  averaging 
of  the  phase  differenees. 
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Strehl  I  Strehl 


Strehl  vs.  Rytov  and  Incoherent  Beacon  Size  for  Least  Squares  Reconstruction 
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5.2.2  Coherent  Extended  Beacon  Performance 

The  coherent  extended  beacons  were  created  using  arrays  of  independent  point 
sources.  Figure  25  shows  the  beacon  intensity  profiles  at  the  receiver  for  five  coherent 
beacon  diameters.  These  profiles  are  for  a  50  km  vacuum  propagation  and  give  some 
insight  into  the  quality  of  their  Gaussian  approximation.  Although  created  with  arrays  of 
point  sources,  the  profiles  show  good  Gaussian  shape  and  the  tail  values  are  very  close  to 
zero,  as  desired.  The  intensity  values  for  each  beacon  have  been  normalized  in  the  plot  to 
provide  for  easy  comparison.  Since  the  reconstructors  rely  on  phase  information  only, 
beacon  intensity  normalization  is  not  required  during  simulation. 


Figure  25:  Normalized  Intensity  vs.  Pixels  for  5  Coherent  Extended  Beacons 


Figure  26  shows  the  performance  of  Goldstein  and  least  squares  algorithms  using 
the  5cm  coherent  beacon  constructed  from  a  point  source  grid.  For  reference,  this  figure 
also  includes  the  Strehl  ratios  obtained  using  an  ideal  point  source  beacon  during  the 
wave-front  sensor  study  in  Chapter  4.  The  Goldstein  point  source  data  line  is  from 
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Figure  7  and  we  see  its  wraparound  problem  at  the  stronger  turbulenee  strengths  sinee  the 
atmospherie  grid  was  not  inereased  from  256x256.  By  eomparison,  the  problem  has  been 
removed  from  the  extended  beacon  data  through  the  use  of  the  512x512  atmospheric  grid 
as  seen  for  the  5cm  beacon.  Similarly,  the  upturn  in  least  squares  performance  at  C„  = 
5x10'  m'  has  been  removed  from  the  extended  beacon  data.  A  slight  amount  of 
wraparound  is  present  for  Goldstein’s  algorithm  at  the  last  two  data  points  of  the 
extended  beacon  results,  but  it  is  not  severe.  Since  5cm  is  less  than  the  diffraction  limit 
of  our  system,  we  see  that  performance  is  close  to  that  obtained  with  a  point  source. 

Figure  27  shows  the  same  information  as  Figure  26,  except  it  is  plotted  versus 
Rytov  and  does  not  include  the  point  source  information. 

Figure  26  through  Figure  31  show  system  performance  for  the  5,  10  and  15cm 
coherent  beacons  plotted  against  C„  and  Rytov.  Comparisons  of  reconstructor 
performance  versus  beacon  size  are  provided  in  Figure  32  through  Figure  36.  We  see 
that  branch  point  detection  is  still  possible  in  the  case  of  coherent  extended  beacons. 

The  advantage  of  Goldstein’s  algorithm  in  strong  turbulence  is  evidenced  in 
Figure  26  through  Figure  31.  For  each  beacon  size,  Goldstein  outperforms  least  squares. 
The  only  exception  is  at  C„  =  1x10'  m'  for  the  15cm  beacon  as  seen  in  Figure  30. 

The  relative  advantage  of  Goldstein’s  algorithm  is  reduced  with  increasing  beacon  size, 
but  is  still  significant  -  1.5  orders  of  magnitude  for  5  cm  scaling  to  slightly  over  1  order  at 
10cm  and  about  1  order  of  magnitude  at  15cm.  The  total  degradation  due  to  beacon 
extent  also  appears  to  be  fairly  uniform  across  the  turbulence  range,  as  the  overall  shape 
of  the  Strehl  ratios  is  maintained  for  each  beacon  size  and  reconstructor. 


61 


Strehl 


Figure  26:  Strehl  vs.  Cn2  for  Ideal  Point  Source  and  5cm  Coherent  Beacon 


0  12  3 

Rytov 


Figure  27:  Strehl  vs.  Rytov  for  5cm  Coherent  Beacon 
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One  result  of  this  shift  is  that  Goldstein  and  least  squares  perform  considerably 
worse  than  uncompensated  below  C„  =4x10'  m'  with  the  15cm  beacon  (Figure  30). 
This  is  due  to  coherent  speckle,  which  creates  field  perturbations  unrelated  to  the 
atmospheric  turbulence.  The  field  sensor  measures  both  the  atmospheric  and  coherent 
speckle  effects  together.  This  creates  errors  in  both  least  squares  and  Goldstein’s 
reconstructions.  Goldstein’s  can  be  significantly  affected  because  the  coherent  speckle 
also  creates  branch  points  in  the  field,  changing  the  placement  of  branch  cuts.  In  weak 
turbulence,  the  phase  error  cost  for  unwrapping  around  these  coherent  speckle  branch 
cuts  is  high,  and  the  reconstructed  phase  error  dominates  the  performance. 

Figure  32  and  Figure  33  show  Goldstein  performance  for  each  beacon  size.  Once 
again,  the  relatively  uniform  nature  of  the  extended  beacon  size  effect  is  visible.  This  is 
logical,  in  that  the  coherent  speckle  is  related  to  the  beacon  size,  which  is  fixed  for  each 
data  set.  We  also  note  that  some  wraparound  is  present  in  all  three  data  sets  for  the  last 
three  C„  values.  The  wraparound  effect  is  minor  and  appears  proportional  to  the  nominal 
Strehl  ratio  at  each  point.  Re-running  these  9  data  sets  using  a  1024x1024  atmospheric 
grid  could  eliminate  this  issue,  but  this  was  deemed  impractical  due  to  time  constraints. 

For  turbulence  strengths  below  4x  10'^^  the  15cm  beacon  returns 
performance  less  than  uncompensated.  Above  4x10'  m'  ,  we  note  that  the  10cm 
beacon  performance  begins  to  approach  that  of  the  15cm  beacon.  This  could  be 
occurring  for  various  reasons.  One  possibility  is  that,  as  C„  increases,  decreases  and, 
as  the  10  and  15cm  beacons  become  less  resolvable,  they  begin  to  appear  more  similar  to 
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Figure  28:  Strehl  vs.  C„  for  10cm  Coherent  Beacon 


Figure  29:  Strehl  vs.  Rytov  for  10cm  Coherent  Beacon 
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Figure  30:  Strehl  vs.  C„  for  15cm  Coherent  Beacon 


Figure  31:  Strehl  vs.  Rytov  for  15cm  Coherent  Beacon 
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Figure  32:  Strehl  vs.  for  Goldstein  Algorithm 


Figure  33:  Strehl  vs.  Rytov  for  Goldstein  Algorithm 
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the  system.  Another  possibility  is  that,  as  the  turbulenee  strength  inereases,  the 
eorrelation  between  the  eorreetion  and  the  atmosphere  deereases  and  both  beaeons  begin 
to  eonverge  on  the  same  partially  eorrelated  performanee. 

The  order  of  magnitude  differenee  between  Goldstein  performanee  and 
uneompensated  at  IxlO"'^  m'^^^  indieates  that  the  10  and  15em  beaeon  eorreetions  are  not 
eompletely  deeorrelated  to  the  atmosphere.  In  strong  turbulenee,  Goldstein’s 
performanee  will  degrade  to  essentially  uneompensated  performanee  when  spatially 
deeorrelated  with  the  ineident  atmosphere,  as  seen  in  Figure  10-  Figure  15  from  the  time- 
delayed,  baseline  geometry. 

While  least  squares  performanee  shows  similar  degradations  to  Goldstein  below 
=  2x10'  m'  ,  it  appears  dominated  by  the  turbulenee  strength  above  this  value.  This 
is  illustrated  in  Figure  34,  where  the  performanee  of  the  10  and  15em  beaeons  generally 
approximate  that  of  uneompensated  performanee.  Figure  35  plots  the  same  information 
as  Figure  34  against  Rytov.  Like  Goldstein’s,  when  least  squares  eorreeted  performanee 
is  essentially  the  same  as  uneompensated  in  strong  turbulenee,  this  generally  indieates 
that  the  eorreetion  is  spatially  deeorrelated  from  the  atmosphere  at  the  time  of  baek 
propagation.  We  also  note  that  the  15em  beaeon  performs  worse  than  uneompensated 
below  C„  =  4x  10'  m'  ,  as  it  did  with  Goldstein.  Figure  36  provides  a  useful 
eomparison  of  least  squares  to  Goldstein  in  strong  turbulenee.  The  best  least  squares  ease 
at  5em  is  outperformed  by  a  lOem  Goldstein  reeonstruetion  for  Rytov  >0.5  and  is  also 
passed  by  the  15em  Goldstein  reeonstruetion  for  Rytov  >1.1. 
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Figure  34:  Strehl  vs.  C„  for  Least  Squares  Algorithm 
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We  can  also  plot  performance  versus  6^  which  can  give  insight  into  the  affect  of 
anisoplanatism.  Since  our  propagation  distance  was  fixed  at  Z  =  50  km,  the  approximate 
isoplanatic  patch  is  given  by  LOo-  Figure  37  shows  the  Strehl  ratio  versus  L6o  for 
Goldstein’s  reconstruction.  The  performance  versus  LO^  for  least  squares  is  given  in 
Figure  38.  For  these  plots,  turbulence  strength  increase  from  right  to  left.  The 
performance  trends  are  closely  related  to  those  seen  for  Rytov  and  as  6o  decreases, 
performance  degrades. 
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Strehl  Strehl 


Figure  37:  Strehl  vs.  LO^  for  Goldstein's  Algorithm 


Figure  38:  Strehl  vs.  LOo  for  Least  Squares 
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Figure  39  through  Figure  44  show  the  mean  Strehl  ±la  for  eaeh  beaeon  eompared 
against  the  mean  uneompensated  Strehl  ±la.  For  eaeh  plot,  the  reeonstrueted 
performanee  is  shown  using  a  dashed  line  and  uneompensated  using  a  dotted  line.  The 
uneompensated  Strehl  varianee  deereases  with  inereasing  turbulenee  strength.  This  is  a 
reasonable  result,  eonsidering  our  baek-propagated  souree.  We  start  with  a  uniformly 
illuminated  laser  whieh  is  foeused  on  the  beaeon  plane  with  our  transmit/reeeive  aperture. 
In  weak  turbulenee,  this  beam  is  affeeted  primarily  by  random  tilt,  whieh  jitters  the 
foeused  spot  around  the  optieal  axis  in  the  beaeon  plane  -  eausing  a  fairly  large  varianee 
over  multiple  independent  samples.  Once  strong  scintillation  begins  to  occur,  the  central 
spot  is  broken  up  into  smaller  regions,  or  is  destroyed  altogether  by  the  turbulence.  This 
results  initially  in  the  energy  being  spread  more  uniformly  at  the  beacon  plane  with  much 
smaller  peak  intensities  in  the  “hot  spots”,  and  hence  yields  a  low  mean  Strehl  and  low 
variance. 

For  a  Goldstein-corrected  back-propagated  beacon,  we  see  a  different  behavior. 
The  Strehl  variance  tends  to  stay  nearly  constant  and  actually  increase  slightly  for  large 
Cn .  According  to  a  1998  Yura  and  Fried  paper  [34],  the  Strehl  variance  is  dominated  by 
the  Rytov  variance  for  well-corrected  systems.  Since  Goldstein’s  algorithm  is  able  to 
maintain  “good”  correction  in  strong  turbulence,  we  see  the  upward  trends  and  large 
variance  in  the  Strehl  performance,  especially  for  the  5cm  beacon.  It  appears  that  the 
variance  is  reduced  as  the  mean  Strehl  performance  decreases,  as  shown  in  Figure  40  and 
Figure  41. 
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Figure  40:  Strehl  ±la  for  Goldstein  lOem 


72 


Figure  41:  Strehl  ±la  for  Goldstein  15cm 


Figure  42-Figure  44  show  the  least  squares  mean  Strehl  ±la.  In  strong 
turbulence,  the  least  squares  Strehl  variance  trends  with  the  uncompensated  variance. 

This  is  best  illustrated  in  Figure  39,  where  a  clear  break  is  present  between  C„  =  1x10' 
and  2x10'  m'  .  Looking  back  to  Figure  34,  this  region  is  also  where  the  least  squares 
performance  tends  to  “break”  towards  uncorrected.  Similar  trends  are  present  in  for  the 
10cm  and  15  cm  beacons,  but  initial  variances  are  lower  due  to  the  decreased  mean  Strehl 
performance.  Based  on  the  strong  overlap  between  corrected  and  uncompensated 
performance  for  the  15cm  beacon,  an  argument  could  be  made  that  least  squares  was 
returning  essentially  uncorrelated  corrections  for  the  entire  turbulence  range.  Conversely, 
Goldstein  mean  Strehl  minus  la  is  better  than  uncompensated  +  la  for  all  but  the 
strongest  turbulence  value  on  the  15cm  beacon. 
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Figure  43:  Strehl  ±la  for  Least  Squares  lOem 
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Figure  44;  Strehl  ±la  for  Least  Squares  15cm 


This  chapter  has  presented  the  simulation  results  for  this  thesis.  We  have  seen 
that  branch  point  algorithms  are  appropriate  for  extended  coherent  beacons  and  point 
source  beacons,  but  do  not  function  with  incoherent  beacons.  The  baseline  studies 
provide  insights  into  the  time-delay  and  noise  performance  of  the  various  reconstructors. 
Chapter  6  will  present  final  conclusions  and  recommendations  for  future  research. 
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VI.  Conclusions 


While  incoherent  extended  beacons  appear  incompatible  with  branch  point 
algorithms,  this  research  demonstrates  the  potential  of  branch  point  tolerant 
reconstructors  to  outperform  traditional  least  squares  techniques  with  point  source  or 
coherent  extended  beacons.  Although  both  least  squares  and  Goldstein  reconstructor 
types  suffer  performance  degradations  due  to  an  extended  coherent  beacon,  Goldstein’s 
algorithm  still  maintains  its  ability  to  deal  with  branch  points  and  calculates  much  better 
phase  corrections  than  least  squares.  It  is  also  clear  that  the  limiting  beacon  size  for  these 
algorithms  is  around  twice  the  diffraction  limit  of  the  optical  system. 

6.1  Importance  of  Modeling  Considerations 

The  importance  of  proper  modeling  discussed  in  Chapter  4  cannot  be  understated. 
Nearly  half  of  the  total  man-hours  involved  in  this  research  focused  on  building  and 
debugging  models  -  particularly  for  the  incoherent  extended  beacon.  Modeling  issues 
that  return  “good”  results  for  certain  scenarios,  but  not  for  others,  can  be  particularly 
perplexing.  This  was  especially  the  case  while  attempting  to  construct  the  speckle 
average  incoherent  beacon  model  using  average  phase.  The  model  appeared 
exceptionally  good  in  vacuum,  but  fell  apart  under  turbulence  propagations.  Finding  the 
answer  required  gaining  a  better  understanding  of  what  was  really  happening  in  the 
simulation  environment.  The  shift  to  average  phase  differences  provided  the  proper 
solution. 

Setting  proper  simulation  parameter  values  begins  with  “rules-of-thumb” 
provided  by  the  simulation  authors,  but  must  be  backed  up  with  experimental  data,  such 
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as  the  required  move  to  5 12x512  grids  in  the  larger  turbulence  when  the  rule-of-thumb 
equations  predicted  256x256  would  be  sufficient.  The  comparison  of  theoretical  Rytov 
variance  to  the  experimental  log  amplitude  variance  in  Figure  5  also  provided  a  good 
indication  that  the  atmospheric  modeling  was  adequate. 

It  was  also  critical  to  establish  a  performance  baseline  by  recreating  the  results  of 
Roggemann  and  Koivunen  [4]  under  ideal  assumptions.  This  validated  the  reconstructor 
modeling  in  combination  with  our  turbulence  model  and  provided  the  basis  for 
performance  comparisons  under  non-ideal  assumptions. 

6.2  Baseline  Geometry  Studies 

The  baseline  geometry  provided  important  insights  into  reconstructor  behavior 
under  various  conditions.  It  confirmed  that  Goldstein’s  algorithm  provides  a  good  upper- 
bound  on  branch  point  tolerant  reconstructor  performance  for  both  zero  and  non-zero 
correction  time  delays.  An  examination  of  least  squares  performance  vs.  Rytov  also 
gives  a  good  sense  of  where  branch  points  begin  to  dominate  turbulence  effects.  The 
noise  study  demonstrated  that  branch  point  cutoffs  could  be  modified  to  accommodate 
measurement  error  -  but  that  the  choice  of  cutoff  becomes  increasingly  difficult  with 
increasing  noise  strength. 

Perhaps  the  most  significant  results  from  the  baseline  geometry  are  the  Strehl  vs. 
field  sensor  resolution  measurements.  Figure  6  provides  a  good  bound  on  branch  point 
reconstructor  performance  vs.  field  sensor  dimensions.  Certainly,  a  worthwhile  practical 
implementation  of  a  branch  point  reconstructor  will  require  at  least  a  16x16  field  sensor 
for  proper  branch  point  detection  in  strong  turbulence.  A  32x32  resolution  would  be 
even  better. 
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6.3  Incoherent  Extended  Beacon  Study 

Incoherent  extended  beaeons  and  branch  point  algorithms  do  not  mix.  As  shown 
in  Figure  24,  beacon  incoherence  frustrates  branch  point  detection.  Since  the  braneh 
point  location  information  is  critical  to  proper  phase  retrieval,  the  hidden  phase  algorithm 
reverted  to  a  least  squares  algorithm.  The  loss  of  branch  point  information  is  most  likely 
due  to  the  phase  averaging  conducted  as  part  of  the  incoherent  assumption. 

As  diseussed  in  Section  3,  incoherence  is  generally  assumed  for  illuminated 
targets.  For  least  squares  algorithms,  incoherenee  is  desirable  to  avoid  the  effeets  of 
eoherent  speekle,  which  induce  field  perturbations  unrelated  to  atmospherie  turbulence. 
This  is  demonstrated  in  the  data  where  the  performanee  degradation  due  to  a  16cm 
incoherent  beaeon  is  less  severe  than  the  16cm  coherent  beacon  for  least  squares. 

6.4  Performance  Bounds  with  Coherent  Extended  Beacons 

The  5cm  coherent  beacon  case  demonstrated  that  ideal  performanee  ean  be 
achieved  with  finite  beacon  extent  if  that  extent  is  unresolvable  to  the  system.  As  beacon 
extent  exceeds  the  diffraetion  limit,  Strehl  performanee  begins  to  suffer  immediately. 

One  interesting  observation  is  the  near  uniformity  of  the  beacon  degradation  effeet  across 
the  turbulenee  range,  as  best  seen  in  Figure  33  and  Figure  35.  Because  of  this  uniform 
degradation,  there  is  motivation  to  stop  AO  corrections  in  weak  turbulence  if  the  beacon 
grows  beyond  a  certain  size.  The  erossover  point  between  correeted  and  uncompensated 
performance  in  terms  of  C„  is  nearly  the  same  for  both  least  squares  and  Goldstein’s 
algorithm. 

The  uniform  degradation  effect  is  dominated  by  coherent  speckle.  While  beacon 
coherence  appears  to  be  required  for  true  braneh  point  detection,  extended  eoherent 


78 


beacons  also  induce  branch  points  which  are  not  related  to  the  atmospheric  turbulence. 
Since  the  coherent  speckle  is  directly  related  to  the  beacon  size  and  not  the  turbulence 
strength,  the  effect  is  somewhat  uniform  vs.  C„  or  Rytov.  Where  few  true  branch  points 
exist,  the  non-atmospheric  field  perturbations  cause  large  phase  reconstruction  errors  . 
This  penalty  is  less  significant  for  strong  turbulence  where  the  benefits  of  detecting  the 
turbulence  branch  points  outweigh  the  performance  cost  of  also  correcting  for  the 
coherent  speckle  branch  points. 

The  ideal  performance  bounds  of  branch  point  algorithms  provide  a  strong 
incentive  for  implementation,  which  is  further  strengthened  after  examining  extended 
coherent  beacon  performance.  Goldstein’s  algorithm  was  able  to  outperform  diffraction 
limited  least  squares  above  0.5  Rytov  with  a  10cm  beacon,  and  above  1.1  Rytov  with  a 
15cm  beacon,  as  shown  in  Figure  36.  This  results  in  a  nearly  a  full  order  of  magnitude 
performance  enhancement  at  the  strongest  turbulence  strength. 

The  Strehl  variances  also  gave  a  good  indication  of  the  performance  of  each 
reconstructor.  Goldstein  Strehl  variances  showed  dependence  on  Rytov  number,  a 
relationship  mathematically  derived  by  Yura  and  Fried  [34]  for  well-corrected  systems. 
Least  squares  only  showed  this  dependence  at  the  lower  turbulence  strengths. 

6.5  Impact  to  Airborne  Laser  and  Other  Directed  Energy  Systems 

This  thesis  has  direct  implications  for  future  development  on  the  Airborne  Laser 
or  derivative  systems.  Improvements  in  reconstructor  performance  offer  to  increase  the 
range  of  ABL  and  reduce  the  fuel  requirements  for  the  high-energy  laser.  This  increase 
in  the  engagement  envelope  has  obvious  dividends  for  force  protection. 
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For  any  choice  of  ABL  reconstructor,  smaller  beaeons  are  always  preferred  to 
maximize  performanee.  The  major  differenees  between  reeonstruetors  depend  on  the 
type  of  beaeons  encountered.  For  point  souree  beaeons,  the  ehoiee  of  reeonstruetor 
depends  primarily  on  whether  braneh  points  are  expeeted  or  not.  In  weak  turbulenee,  the 
use  of  least  squares  is  entirely  appropriate,  but  the  ehoiee  of  a  braneh  point  algorithm  is 
elearly  motivated  in  strong  turbulenee.  In  the  ease  of  extended  beaeons,  the  beaeon 
eoherenee  plays  a  large  role  in  the  ehoiee  of  reeonstruetor. 

While  the  ABL  beaeon  is  generally  assumed  to  be  ineoherent  for  modeling 
purposes,  AFRL/DE  has  found  that,  in  fact,  it  is  better  modeled  as  a  partially  eoherent 
beaeon.  This  partial  eoherenee  eould  provide  the  neeessary  braneh  point  information  to 
implement  a  branch  point  reeonstruetor  on  the  system.  While  the  resulting  performanee 
would  be  less  than  the  limiting  ease  of  full  eoherenee,  it  eould  well  outperform  least 
squares  reeonstruetion  in  strong  turbulenee  and  is  worthy  of  investigation. 

Future  systems  eould  also  benefit  from  enhaneed  reeonstruetor  performanee  in 
many  ways.  Praetieal  engagement  range  eould  be  achieved  with  a  lower  power  kill  laser. 
Robustness  to  an  extended  beaeon  eould  reduee  the  eomplexity  of  the  illumination 
system  and  henee  its  eost  or  weight.  Together,  these  improvements  could  reduce  the 
paekage  size  of  an  AO  system  and  make  it  more  praetieal  for  taetieal  deployment  on  non- 
eargo  airframes. 

Future  systems  eould  be  designed  to  maximize  beaeon  eoherenee  for  operations  in 
strong  extended  turbulence.  This  would  provide  the  maximum  benefit  to  a  braneh  point 
reeonstruetor.  Weak  turbulenee  performanee  eould  be  kept  elose  to  point  souree 
performanee  as  long  as  the  beacon  size  was  well  eontrolled. 
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6.6  Future  Research 


As  with  any  good  scientific  endeavor,  the  attempt  to  answer  a  few  questions 
opens  the  door  to  many  more.  Follow-on  topies  can  be  roughly  divided  into  three  groups. 
The  first  are  topics  related  to  the  coherent  extended  beaeon,  modeled  using  the  point 
source  array.  Second  are  topics  involving  general  aspects  of  AO  modeling  encountered 
in  this  work,  and  third,  topics  which  push  beyond  the  coherent  extended  beacon  into 
larger  non-ideal  assumptions. 

6.6.1  Possible  Topics  on  Gaussian  Extended  Sources 

Application  of  other  branch  point  tolerant  reconstructors  to  the  extended  beacons 
would  be  very  interesting.  Particularly,  application  of  Tyler’s  slope  discrepancy  [5]  or 
Roggemann’s  intensity  weighted  least  squares  algorithm  [17]  would  provide  performance 
estimates  for  realizable  reconstructor  algorithms.  These  techniques  do  not  depend  on 
branch  point  locations  -  so  perhaps  they  would  perform  better  with  incoherent  beacons. 
The  main  difficulty  in  such  a  project  would  be  implementation  of  the  reconstructor  and 
choice  of  field  sensor.  The  beacon  and  atmospheric  models  could  likely  remain 
unchanged.  The  use  of  a  practical  wave-front  sensor  would  also  allow  for  a  noise  study 
for  any  of  these  cases. 

Fried  has  also  developed  an  exponential  class  reconstructor,  which  constructs 
wrapped  phase  from  phase  differences  [35].  This  reconstructor  has  been  formulated  to 
maintain  the  branch  point  information  while  exhibiting  some  degree  of  white-noise 
rejection.  This  reconstructor  could  provide  an  important  front  end  to  Goldstein’s 
algorithm  and  allow  it  to  be  implemented  in  practice. 
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It  will  also  be  important  to  examine  reconstructor  performance  with  partially 
coherent  beacons.  Beacon  modeling  will  be  a  large  portion  of  such  a  study.  Since  most 
practical  incoherent  illuminator  systems  actually  demonstrate  some  degree  of  partial 
coherence,  branch  point  reconstruction  algorithms  may  still  be  applicable.  As  discussed 
in  Section  6.5,  the  Airborne  Laser  is  one  such  system. 

6.6.2  General  Extended  Beacon  Topics 

The  next  major  step  in  extended  beacon  modeling  is  to  build  a  model  for  an 
extended,  scintillated  beacon.  This  type  of  beacon  occurs  for  cases  where  the  target  is 
illuminated  by  a  beacon,  which  first  propagates  through  extended  turbulence  to  the  target. 
Wave-front  sensing  on  such  a  beacon  is  expected  to  be  very  difficult,  and  new  techniques 
will  likely  be  required  to  achieve  acceptable  performance. 

One  of  the  most  likely  techniques  will  be  wave-front  sensor  compensation  based 
on  beacon  intensity  estimation.  Estimating  the  intensity  of  an  extended,  scintillated 
beacon  will  likely  require  some  type  of  blind  deconvolution  algorithm.  The  problem  is 
made  particularly  difficult  because  the  problem  is  doubly  stochastic  -  a  random 
fluctuating  beacon  propagating  through  a  random  fluctuating  atmosphere.  Most  blind 
deconvolution  algorithms  assume  a  constant  “target”  or  beacon  and  iterate  over  multiple 
“looks”  [36]  which  is  not  appropriate  for  this  problem. 

Assuming  a  workable  beacon  estimator  is  developed,  the  next  question  is,  how  to 
use  this  information.  Clearly,  extended  beacons  corrupt  wave-front  sensor  measurement. 
The  challenge  is  to  derive  a  relationship  for  this  corruption  that  could  be  used  to  pre¬ 
correct  wave-front  sensor  slope  measurements  prior  to  phase  reconstruction. 
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Alternatively,  perhaps  a  phase  retrieval  algorithm  could  be  formulated  to  account  for  the 
beacon  intensity  estimate  in  some  way. 

6.6.3  Topics  Derived  from  General  AO  System  Observations 

One  research  topic  could  be  a  parameter  study  of  the  atmospheric  coherence 
diameter,  isoplanatic  angle,  Rytov  variance,  and  integrated  C„  (ro ,  6o,  andM^)  vs. 

Strehl  variance  to  prove  Yura  and  Fried’s  relationship  [34]  for  Strehl  variance 
experimentally.  This  could  be  done  using  a  point  source  beacon  with  the  Goldstein  and 
least  squares  reconstructors  over  a  wider  range  of  turbulence  strengths.  Another  avenue 
would  be  to  probe  the  relationship  between  mean  Strehl  and  variance  and  the  spatial 
correlation  of  the  correction  vs.  atmosphere.  This  could  be  done  using  WaveTrain’s  wind 
models  and  time  delay  features. 

During  the  course  of  working  with  Goldstein’s  algorithm,  it  was  discovered  that 
disabling  branch  point  finding  only  reduced  point  source  beacon  performance  about  5%. 
This  small  performance  degradation  is  not  easily  explained,  especially  for  strong 
turbulence,  but  was  obtained  for  noise-free  point  source  scenarios.  In  this  case,  it  may  be 
that  the  SNR  is  sufficiently  high  to  render  any  path  dependent  unwrapping  superior  to 
least  squares.  An  obvious  question  is  to  determine  if  the  SNR  is  truly  the  dominant  effect 
and  quantify  its  relationship  to  the  unwrapping  performance. 

6.7  Final  Thoughts 

Improvements  in  adaptive  optics  continue  to  create  new  applications  and 
opportunities.  While  the  directed  energy  and  laser  communication  applications  are  the 
current  drivers  behind  extended  beacon  and  turbulence  research,  future  applications  will 
likely  involve  much  smaller  scales.  As  practical  fabrication  of  MEMs  devices  increases. 
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there  will  be  a  push  to  develop  miniaturized  AO  systems.  There  is  also  the  desire  to  use 
less  expensive  components  in  laser  projection  systems,  wave-front  sensors,  and  mirrors. 
Reducing  the  dependence  on  beacon  extent  will  reduce  requirements  on  illuminator  lasers 
and  optics  -  helping  to  lower  the  cost  of  AO  systems. 
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